M Liévin
email: lievin@tirf.inpg.fr

F Luthon
email: luthon@tirf.inpg.fr

UNSUPERVISED LIP SEGMENTATION UNDER NATURAL CONDITIONS

Unsupervised lip segmentation under natural conditions

INTRODUCTION

It is commonly observed that visual information provides a precious help to the listener under degraded acoustical conditions [START_REF] Benoît | A set of french visemes for visual speech synthesis[END_REF]. The motivation of the present work is to extract visual information for automatic speech recognition (ASR), videoconferencing and speaker's face synthesis under natural lighting conditions with few assumptions. Some approaches proposed in this area are based on gray level analysis (e.g. Luettin in [START_REF] Stork | Speechreading by Humans and Machines[END_REF]). Others use color analysis but need to determine optimal values of some parameters (e.g. Coianiz in [START_REF] Stork | Speechreading by Humans and Machines[END_REF]). Strong assumptions are required on the skin hue parameters and the mouth location [START_REF] Wark | A synthetic approach to automatic lip feature extraction for speaker identification[END_REF], therefore the skin hue region is often determined manually beforehand.

The previous work [START_REF] Liévin | Lip features automatic extraction[END_REF] used a segmentation to locate the mouth before estimating lip geometrical features, some of the segmentation parameters were determined beforehand.

Here, an algorithm is proposed for unsupervised lip shape extraction and mouth location under natural conditions, the requirement being that a micro-camera is mounted on a light helmet worn by the speaker so that it is fixed w.r.t. the head.

The RGB video sequence (8 bits/color/pixel) contains the region of the face spanning from chin to nostrils. The purpose of the process is to obtain the mouth shape using red hue label fields and motion information. The processing is divided into three stages:

1: Logarithmic color-space transform, RGB to HI. 2: of sequence dependant parameters:

Computation of the mean value of the lip hue H lip .

Estimation of the noise on the motion information. 3: Spatiotemporal segmentation of the lip and ROI estimation.

PARAMETER ESTIMATION

Logarithmic color transform

Color-based approaches often use color angles methods (HSI) for illuminant-invariant recognition. Color shifts can be well categorized with angles if camera sensors are sufficiently narrowband. But, in our application, we deal with a mono-CCD camera which gives poor results with angular transforms (noisy conditions). Moreover, G and B channels seem to be correlated in the red region (where R is preponderant). From the RGB color space, we use only two dimensions R and G under the assumption that red prevails in face areas and specially in lip areas. We define the chro- To obtain a robust hue observation to the lighting conditions, we compute the hue in a mathematical framework based on a logarithmic image processing model [START_REF] Jourlin | Image dynamic range enhancement and stabilization in the context of the logarithmic image processing model[END_REF]. The intensity I of an image is represented by its associated gray tone function i (Eq. 1). This model satisfies the saturation characteristics of the human visual system and is justified from a physical point of view. Specific algebraic and functional operations are redefined in a vectorial structure, like and as respectively the addition and the opposite of a gray tone function. The difference between g and f, respectively logarithmic tone of the intensity G and F ,isgivenEq.

2. i = M 1 , I I 0 (1) f g = M f , g M , g (2) h = g r = M 1 , G R (3)
Considering the illuminance I 0 close to the maximum value of white M , the logarithmic transform becomes i = M ,I.

We define h as the logarithmic hue tone of H, difference of g and r, logarithmic color tone of G and R (Eq. 3). The logarithmic difference becomes a ratio between R and G components H = M G R . Finally, from the RGB color space, a HI logarithmic color space (Fig. 2) is defined considering M =256and the intensity I as the mean value of the R,G and B components (Eq. 4). the corresponding hue sequence.

H =256 G R and I = R + G + B 3 (4)

Observations

To detect lip regions, motion information is combined with red hue. From the HI color space, two kinds of observations o are derived, defined to be in the same range 0 255 as the image quantification (8 bits). First, a hue observation hs consists in filtering the hue value Hs at pixel s with a parabola centered on the mean value of lip hue H lip with a standard deviation of the hue value H (Eq. 5). H 16 (5) The notation 1 condition denotes a binary function which takes the value 1 is the condition is true, 0 otherwise.

Second, a temporal observation fds is defined as unsigned difference between the luminance of two consecutive images (Eq. 6). Is represents the intensity (or luminance) at pixel s. fds=jI t s , I t,1 sj (6)

Hue and motion estimation

The hue segmentation needs three estimated parameters to be unsupervised : H lip , H , h . Due to the chosen expression of hs, the tresholded hue field is defined by Eq. 7, expressing the link between H and h . hs h ,jHs , H lip j H p 256 , h (7)

The hue histogram fp i;h ;i 2 0 255 g is an useful representation of the hue distribution over the image. We can detect two modes, one for the face, the second for the lip. But, in natural conditions (no make-up), the two modes are mixed (Fig. 3). In order to estimate H lip accurately, a spe- cific hue H h f ace;lip is defined in Eq. 8. H lip is defined by the Eq. 9.

H h S = cardS h X i2 h p i;h S i (8)

H lip = H h lip S s (9)
h corresponds to the appropriate interval and S s represents the image S after face segmentation. The processing respects the following steps:

1: Estimate H face using fp i;h ;i 2 0 255 g computed over the hue image. 2: Segment the first hue image with a basic spatial MRF segmentation via the hue transform observation The algorithm requires an appropriate threshold fd to suppress the camera noise without cutting significant temporal changes. In the previous work [START_REF] Liévin | Lip features automatic extraction[END_REF], this threshold was determined before segmentation by hand. We compute here the entropy E fd S over an image difference (Eq. 11). This gives the level of noise from which we can deduce the value of fd (Eq. 12). The thresholded motion field is then defined by fd fd . The thresholded fields appear non homogeneous and noisy (Fig. 4). Therefore, we need a statistical relaxation to segment more accurately the lip.

THE SEGMENTATION ALGORITHM

The spatiotemporal MRF framework

From these two thresholded observations, four initial labels (a 0 ,a 1 ,b 0 ,b 1 ) are derived, for coding four pixel classes: pixels with ( 1 ) (resp. without ( 0 )) motion, belonging (a) (resp. not belonging (b)) to red hue areas. This label field is supposed to follow the main MRF (Markov Random Field) property related to a spatiotemporal neighborhood structure, i.e. the label l s of the current pixell s depends only on the labels of its spatiotemporal neighbors n.

Maximizing the A Posteriori probability (MAP criterion) of the label field is equivalent to minimizing a global energy function [START_REF] Geman | Stochastic relaxation, gibbs distributions and the bayesian restoration of images[END_REF]:

W S= X o2ffd;hg U o S+ :U m S (13)
where U o and U m represent respectively the attachment energies (expressing the link between labels and observations, Eq. 14) and the model energy (corresponding to spatial and temporal a priori constraints) (Eq. 15) over the image S, is a weighting coefficient between the two energies.

U o S= X s2S o s , o l s 2 2 2 o ( 14
)
where o is an attachment function, mean value of the observation o over S and 2 o is the corresponding variance.

Both are estimated on line.

The apri orimodel energy is defined as a sum of interaction potential functions over the neighborhood:

U m S= X s2S h X n2s V st l n ;l s i (15)
The spatiotemporal potential function V st is defined as the inverse of the Euclidian distance between two neighbors.

The distance integrates two elementary potentials s and t as scale factors (Eq. 16).

V st l n ;l s = s l n ;l s t l n ;l s q t l n ;l s 2 , where ,,,! s; n= x ; y ; t and 2f,1; 0; 1g

The elementary potentials s and t are defined to con- strain the model respectively to spatial homogeneity of labels and temporal homogeneity of hue when no motion is detected (details in [START_REF] Liévin | Lip features automatic extraction[END_REF]).

The relaxation algorithm

The iterative deterministic algorithm ICM (Iterated Conditional Modes) is implemented to compute the minimum energy at each site (Eq. 13 with typ. =20), starting from the thresholded fields as initial label configuration. After a few iterations on the field (less than 10 to respect the stopping criterion for convergence WS=WS 0:05 %), convergence is achieved. One obtains homogeneous red hue and lip motion fields (Fig. 5).

ROI estimation

From lip red hue relevant labels, the ROI is evaluated on line by maximising a cost function ,S on each image (Eq. 10 in [START_REF] Liévin | Lip features automatic extraction[END_REF]) after each step of the relaxation. One each image, the last estimated ROI is increased with a scale factor and used to initialize the current one. The ROI estimation reduces the relaxation time by surrounding the mouth precisely. Moreover, it increases the accuracy of parameter's estimation.

LIP SHAPE EXTRACTION

Different sequences have been tested, some with natural red make-up (To p in Fig. 6), others with poor lighting conditions without any make-up (Bottom in Fig. 6). These results show the robustness of the unsupervised algorithm to the variability of natural conditions. The unsupervised parameter estimation method gives a one pixel mean difference with ground truth measures for the vertical height and the horizontal width of the internal lip opening (Fig. 7). The external shape is unfortunately more elusive but accurate enough to initiate a simple deformable geometrical model.

CONCLUSION

An unsupervised lip segmentation have been successfully applied to several sequences in natural conditions (natural Bunsup: as unsupervised parameter estimation measure; B g:truth as ground truth measure)

images of speaker's face without any particular make-up or lighting). First, the choice of a logarithmic transformation close to the characteristic of the human visual system enables the algorithm to estimate accurately the mean value of lip hue H lip (the speaker's dependant parameter). This transformation is combined with a noise estimation on the frame difference. Second, the spatiotemporal algorithm integrates hue with motion information, improving the quality of contours often elusive on speaker's lips. Finally, the quality of the segmented fields is similar to those obtained with parameters determined beforehand manually [START_REF] Liévin | Lip features automatic extraction[END_REF]. We need to process more sequences to test the robustness of the parameter estimation with more difficult cases, like faces with beard or colored people faces. The proposed algorithm requires less than 10 iterations until convergence (about 2 sec. on a SunUltra1).

Figure 1 :

 1 Figure 1: Left:3 DRGB histogram of a face under natural illumination; Right:2 DRG chromacity histogram. macity histogram (R,G) as the non-normalized projection of the RGB color space. The typical histogram of a face sample is shown Fig. 1. Two regions with a specific angular direction appear.To obtain a robust hue observation to the lighting conditions, we compute the hue in a mathematical framework

Figure 2 :

 2 Figure 2: Top: 5 typical images of luminance sequence; Bottom:

Figure 3 :

 3 Figure 3: Left: histogram of hue image;Right: unsupervised segmentation of hue face.

  p i;o S represents the probability of level i in the observation o over image S.

Figure 4 :

 4 Figure 4: From top to bottom: sequence of temporal observation fd; sequence of red predominance observation h with unsupervised parameter estimation (H lip = 130 ; H = 6); sequence of temporal observation thresholded with unsupervised parameter estimation ( fd =9); sequence of red prevailing observation thresholded ( h =192).

Figure 5 :

 5 Figure 5: From top to bottom: initial labels; label fields after relaxation, The 4 labels are shown in grey levels (from white to black: b1, a1, b0, a0); red hue relevant label images (a0 and a1) superposed with luminance.

Figure 6 :

 6 Figure 6: Top: Sequence of final red hue fields with ROI superposed on the luminance with soft red make-up; Bottom: Sequence of final red hue fields with ROI superposed on the luminance with no lighting supply and no make-up.

Figure 7 :

 7 Figure 7: Internal lip measurement on the sequence partially shown Fig. 5 (B manual as manual parameter setting measure;

3 :

 3 Evaluate H lip using fp i;h ;i 2 0 255 g computed over the segmented hue image. lip and h face, camera dependant parameters, are independant from the speaker and the lighting conditions. They can be estimated by camera calibration. Currently, the range of h is the result of the statistical distribution of manually estimate over caracteristic natural conditions. The selected range for h face;lip is 100 200 and for h lip, 100 150 . This range corresponds to the red predominant region. The equation 7 is then respected with h =192and H =6.

h


Marc Liévin
email: lievin@lis-viallet.inpg.fr

Patrice Delmas
email: delmas@lis-viallet.inpg.fr

Pierre-Yves Y Coulon

Franck Luthon

Vincent Fristot

Automatic lip tracking: Bayesian segmentation and active contours in a cooperative scheme

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

It is commonly observed that visual information provides a precious help to the listener under degraded acoustical conditions [START_REF] Benoît | A set of french visemes for visual speech synthesis[END_REF]. The motivation of the present work is to extract visual information for automatic speech recognition (ASR), videoconferencing and speaker's face synthesis under natural lighting conditions with few assumptions. Some approaches proposed in this area are based on grey level analysis (e.g. Luettin in [START_REF] Stork | Speechreading by Humans and Machines[END_REF]). Others use color analysis but need to determine optimal values of some parameters (e.g. Coianiz in [START_REF] Stork | Speechreading by Humans and Machines[END_REF]). A wide range of papers describe the applications of active contours for lip boundary detection but often focus on inner (e.g. Petajan in [START_REF] Stork | Speechreading by Humans and Machines[END_REF]) or outer lip contours only, rarely both.

Here, an algorithm is proposed for inner and outer lip contour tracking under natural conditions, the requirement being that a micro-camera is mounted on a light helmet worn by the speaker so that it is fixed w.r.t. the head. The RGB video sequence (8 bits/color/pixel) contains the region of the face spanning from chin to nostrils. The purpose of the process is to obtain accurate inner and outer lip borders even if the mouth is closed. A Bayesian segmentation [START_REF] Liévin | Lip features automatic extraction[END_REF] is used as an initialisation step for a snake convergence. The cooperative scheme is divided into three stages:

Preprocessing stage: Logarithmic color-space transform, RGB to HI. Unsupervised spatiotemporal segmentation of lip area and mouth location estimation.

Gateway: Mouth characteristics extraction:

Boundaries and semantic information from the segmented lips. Mouth corners detection.

Automatic Snake processing: Snake initialization from preprocessing results. Convergence of automatic snake (outer) and balloon snake (inner). This work is part of the Labiophone project (ELESA Federation n. 8 (CNRS-INPG)), an advanced audio-visual communication tool (Fig. 1). This project aims at providing a very low bit rate coding communication system, integrating both audio and visual features.

Lip segmentation

Logarithmic color transform

Face features detection is often illuminance dependent. To gain independence from lighting conditions, we compute here a color transform for illuminant-invariant recognition. Angular transforms give poor results in noisy conditions (mono-CCD camera). Therefore, a logarithmic hue transform is defined using G and B channels from the RGB color space.

We compute the hue in a mathematical framework based on a logarithmic image processing model [START_REF] Jourlin | Image dynamic range enhancement and stabilization in the context of the logarithmic image processing model[END_REF]. The intensity I of an image is represented by its associated gray tone

function i = M (1; I I0
). This model satisfies the saturation characteristics of the human visual system and is justified from a physical point of view. Specific algebraic and functional operations are redefined in a vectorial structure. The difference between logarithmic tone of the channels G and R corresponds to the logarithmic hue tone h. With few assumptions (I 0 close to the maximum value of white M ), the logarithmic difference becomes a ratio between G and R components. Finally, from the RGB color space, a HI logarithmic color space is defined (Eq. 1) (Fig. 2).

H =256 G R and I = R + G + B 3 (1)
Figure 2: Top: 5 images of a typical luminance sequence; Bottom: the corresponding hue sequence.

Lip Hue Segmentation

Observations

To detect lip regions, motion information is combined with red hue. From the HI color space, two kinds of observations o are derived (Eq. 2): h(s) consists in filtering the hue value H(s) at pixel s with a parabola centred on the mean value of lip hue H lip with a standard deviation of the hue value H ; fd(s) is defined as the unsigned difference be-

tween
The notation 1 condition denotes a binary function which takes the value 1 if the condition is true, 0 otherwise.

Hue and motion estimation

Hue and motion parameters are estimated automatically. In the previous work [START_REF] Liévin | Lip features automatic extraction[END_REF], these thresholds were determined before segmentation by hand. The hue observation needs three parameters to be estimated: H lip , H , h . For that purpose, the hue histogram is a useful representation of the hue distribution over the image. We can detect two main modes: the first for the skin-lip face, the second for the background. In natural conditions (no make-up), the lip mode and the skin mode overlap (Fig. 3).

In order to estimate H lip accurately, the processing re- spects the following steps:

Estimate H skin from the hue distribution computed over the whole image (Left of the Fig. 3), the only assumption being that the main mode corresponds to the hue skin.

Cluster all pixels respecting the condition given in Eq.

2 with H skin instead of H lip (Middle of the Fig. 3). Evaluate H lip from the hue distribution after discard- ing all pixels belonging to skin mode. The remaining are the lip mode (black of the Fig. 3) and the background mode.

The threshold hue field is then defined by h> h . The algorithm requires an appropriate threshold fd to suppress the camera noise without cutting significant temporal changes. We compute here the entropy E fd (S) over an image. The threshold motion field is then defined by fd > fd with fd (S)=2 E fd (S) .

The thresholded fields appear non homogeneous and noisy. Therefore, we need a statistical relaxation to segment more accurately the lips.

The Segmentation Algorithm

Observations and Labels in an MRF Framework

From these two thresholded observations, four initial labels (a 0 ,a 1 ,b 0 ,b 1 ) are derived for coding four pixel classes: pixels with ( 1 ) (resp. without ( 0 )) motion, belonging (a) (resp. not belonging (b)) to red hue areas. This label field is supposed to follow the main MRF (Markov Random Field) property related to a spatiotemporal neighborhood structure (Fig. 4), i.e. the label l s of the current pixell s depends only on the labels of its spatiotemporal neighbors n. where U o and U m represent respectively the attachment energies (expressing the link between labels and observations, Eq. 4) and the model energy (corresponding to spatial and temporal a priori constraints) (Eq. 5) over the image S, is a weighting coefficient between the two energies.

U o (S)= X s2S o s ; o (l s )] 2 2 2 o (4)
where o is an attachment function, mean value of the ob- servation o over S and 2 o is the corresponding variance.

Both are estimated on line.

The a priori model energy is defined as a sum of interaction potential functions over the neighborhood:

U m (S)= X s2S h X n2 (s) V st (l n l s ) i (5)
The spatiotemporal potential function V st is defined as the inverse of the Euclidian distance between two neighbors.

The distance integrates two elementary potentials s and t as scale factors (Eq. 6).

V st (l n l s )= s (l n l s ) t (l n l s ) q t (l n l s ) 2 ; 2 x +4 2 y + s (l n l s ) 2 2 t [START_REF] Liévin | Lip features automatic extraction[END_REF] where ;;;! (s n)=( x y t ) and 2f;1 0 1g

The elementary potentials s and t are defined to con- strain the model respectively to spatial homogeneity of labels and temporal homogeneity of hue when no motion is detected.

Results and ROI estimation

An iterative deterministic algorithm (ICM : Iterated Conditional Modes) is implemented to compute the minimum energy at each site, starting from the initial label configuration

L 0 t .
From lip red hue relevant labels, the ROI (Region Of Interest) is evaluated on line by maximizing a cost function ;(S) on each image (details in [START_REF] Liévin | Lip features automatic extraction[END_REF]) after each step of the relaxation. The ROI estimation reduces the relaxation time by surrounding the mouth precisely. Moreover, it increases the accuracy of parameter's estimation.

After a few iterations (typ. 10) on the label field, convergence is achieved. One obtains homogeneous red hue and motion lip fields.

Lip red hue labels and ROI

From the final label fields, one can extract lip red hue relevant label (a 0 and a 1 ) (Fig. 6). Those results are shown with a ROI evaluated on line. Several typical sequences have been tested, some with a soft natural red make-up, others with very poor lighting conditions without any make-up. It shows the robustness of the algorithm to the variability of the lighting conditions.

Preprocessing lip contours extraction

As a preprocessing stage, the lip segmentation offers robust information to lighting conditions: automatic mouth location (ROI) (Fig. 6) unsupervised segmentation of the mouth shape semantic information: open/close detection But, the borders seem to appear irregular when tongue or gum are segmented, the segmentation is elusive when close to the mouth corners, the inner contour is not segmented when the mouth is closed. We need then a higher level mouth detection algorithm: active contours. A gateway has been defined to initialise the next stage with relevant data from the preprocessing.

Intermediate processing: a gateway towards active contours

Mouth corners detection

ROI information is used to locate mouth corners and vertical limits of the lips. Areas of darkness occur at the inner border of lips on horizontal mouth transitions (e.g.: upper lips and tooth, tooth and mouth interior, tooth and lower lips). Indeed, the vertical minima of image locates lip frontiers and corners with accuracy.

The mouth corners are estimated with the following steps:

Find the grey level minima pixel over image columns; compute the distribution (with as core (Eq. 7)).

j =4 j (N col ; j) (N col ) 2 (7)
where N col is the number of columns of the image, j the current column, the weighting coefficient j varies from 0 to 1, from the border to the center of the image.

Detect the highest peak and deduce the horizontal symmetry axis of the mouth.

Extract lip corners following the line of minima, from the center of the image to the left (respectively to the right).

Estimate the width and orientation angle of the mouth.

Finally, a good estimate of the width and the orientation of the mouth can be used (Fig. 7).

Lip shape extraction

Starting from the ROI coordinates, a reduced number of edges (e.g. 30) are detected within the inner and the outer borders of the lip shape. These edges are linked together with the mouth corners we found previously to provide an excellent snake initialization (Fig 8 ).

Active contours

Three major problems are classically encountered while using snakes: initialization, parameter estimation and convergence of the algorithm. Initialization is commonly done by hand, close to the object to provide good convergence. Snake evolution is sensitive to parameter values which are usually evaluated manually after several tests. Snake convergence needs a good fitting between its energy and the desired image features.

Xu [START_REF] Xu | Gradien vector flow: A new external force for snakes[END_REF] solves most of the initialization and convergence problems for concave areas by creating a new external force called gradient vector flow. But the diffusion process required too much time and therefore has no quasi real time applications. L.D Cohen [START_REF] Cohen | On active contour models and balloons[END_REF] introduced balloon forces to provide inflation of snakes or push them far away from their initial guess. This approach helps to ensure the snake convergence or gives snake dilatation properties.

The lip tracking algorithm uses non-convex energy minimization to succeed in detecting mouth boundaries. The first energy minimum reached is likely to be a local one. Therefore, a good mouth parameter estimation is an essential step for good results.

Energy minimizing curves

Introduced by Kass and al [START_REF] Kass | Snakes: Actives contours models[END_REF] active contours were designed for interactive interpretation in which the user guides (by external forces modification) the snake near the desired solution. A snake is a parameterised curve v defined (Eq. 8) by its Cartesian coordinates x and y along the curvilinear abscissa s which evolves through the minimization of its functional (Eq. 9).

v(s) = x(s) y(s)] s 0 1] (8) :v(s) ;! Z 1 0 (E int (s)+E ext (s)) ds (9)
The internal energy (Eq. 10) is a second order regularization term derived from Tikhonov ill-posed problems theory. It controls the curve smoothness via weighting parameters and .

controls the snake tension and its curvature. External energy (Eq. 11) represents the fitting of image data to current vector. We focus on lip boundaries. We then decided to use the image gradient to extract edge points. To do so we use a classical gradient filter (such as Sobel or Canny-Deriche). E int (s) = jv 0 (s)j 2 + jv 00 (s)j 2 (10)

E ext (s) = ;jr (G I)(v(s))j 2 (11)
r represents the gradient operator, G the 2D Gaussian kernel and I the current image. This leads us to the classical dynamic scheme (Eq. 12) where I d is the identity matrix, A the Toeplitz snake matrix, V the snake control points vector and 1 the time step coefficient.

V (t)=(A + I d ) ;1 ( V (t ; 1) ; F(V (t ; 1))) (12)

F represents forces derived from external energy.

Higher level information forces such as Distance map or Balloon forces [START_REF] Cohen | On active contour models and balloons[END_REF] are added there.

Adapted snakes

Our automatic snakes integrate and as non spatially constants values. The Toeplitz matrix obtained is not detailed here. Forces calculation is done by bilinear interpolation to reduce numerical instabilities which occur through snake energy minimization. Sampling the snake curve by spline functions during the process enforces a constant distance between snake points. Moreover, it helps moving points trapped by spurious edges. Finally we impose all parameters constant through different images for a given number of snake points.

As mouth corners are finely detected, the snake gets its extremities fixed. That kind of active contours is less unstable than the traditional ones. Convergence is tested after each resampling (classically every 10 iterations). This solves the final oscillation convergence problema of snakes. We authorize a maximum quadratic error (Eq. 13) calculated between two successive sampled snakes. = X i2 0::N;1] jV i (t + t s ) ; V i (t)j (13) where t s represents the sampling time step chosen and N the size of the snake control points vector.

Shape constraint

Our aim is to maintain a mouth shaped snake even without external constraints. Thus, we test non spatially constant snake parameters derived from physical considerations based on mouth geometry. For example, the lower lip contour usually has a curvature that is minimal at the middle and maximal at the corners. We choose a higher coefficient at the middle and a lower one around corners. The same kind of adjustment is applied to the upper lip.

Results

Using the gateway lip shape information to initialize snakes gives us the opportunity to be closer to the desired object. Inner positioning is done in the mouth. Therefore, we should inflate the inner snake to reach inner lips boundaries. Cohen Balloon forces reduce the natural ability of snakes to shrink even without external forces. When the mouth is closed, we simply sample the inner snake along the straight line between lip corners.

We hold all the coefficients constant throughout our tests for both inner and outer snakes. Images from our image database including open and closed mouths from different faces are tested. Outer and inner snakes always reach good boundaries after few iterations, usually less than 100.

The top of Fig. 9 shows five successive frames of a closing mouth. Outer and inner lips are perfectly tracked with a few number of iterations (about 60). The inner snake is capable of detecting small and asymmetric mouth apertures (frames 2 to 4).

The bottom of Fig. 9 points out specific problems of lip tracking. Lips vary in shape from one speaker to another. Benny has thin upper lips. His mouth is longer than the previous speaker. His upper lip is very elusive but our algorithm succeeds in detecting both inner and outer lip frontiers.

Conclusion

An automatic lip contours extraction has been successfully applied to several sequences in natural conditions (natural images of speaker's face without any particular make-up or lighting).

First, the logarithmic transform and the following spatiotemporal Bayesian segmentation dealing with hue and motion information locate and segment the mouth with very few error rate. Then, a preprocessing gateway, combining low level information and mouth corners location, initializes inner and outer active contours. Thanks to the mouth corners location algorithm, the two snakes get their extremities fixed. Whereas snakes were described by their inventors as a semi-automatic process, these preprocessing steps prove that automatic snakes are viable. Finally, with no assumption about the lighting conditions, the lip contours (inner and outer) are extracted with accuracy, even when the mouth is closed or asymmetric. Some problems still occur for a fine detection when tongue or gum are detected. We need also to deal with more difficult cases like colored people or faces with beard (Fig. 10). The proposed algorithm requires less than 3 to 4 seconds per image on a standard 150MHz workstation. Therefore, hardware implementation for both stages, spatiotemporal segmentation and active contours, are currently under study.

Figure 1 :

 1 Figure 1: Context of lipreading: from sequence of speaker images, lip tracking provide parameters for talking face synthesis.

Figure 3 :

 3 Figure 3: From left to right: histogram of hue image (In black: overlap between lip and skin distribution); the corresponding segmentation of skin hue (In black); histogram of hue image when the skin mode is discarded (In black: the lip mode).

Figure 4 :

 4 Figure 4: Left: Spatiotemporal neighbourhood structure with binary cliques c = (s n). s is the current pixel (in black), n is any spatiotemporal neighbour of s (in gray); Right: corresponding elementary cube C xyt Maximizing the A Posteriori probability (MAP criterion) of the label field is equivalent to minimizing a global energy function [3]: W (S)= X o2ffd hg U o (S)+ :U m (S)

Figure 5 :

 5 Figure 5: From top to bottom: sequence of luminance images; initial labels; label fields after relaxation: the 4 labels are shown in gray levels (from white to black: b 1 , a 1 , b 0 , a 0 ); sequence of hue relevant label images (a 0 and a 1 ).

Figure 6 :

 6 Figure 6: Two sequences of final lip hue fields with ROI superposed on the corresponding luminance.

Figure 7 :

 7 Figure 7: From left to right: vertical minima (in white) on a grey level image; the corresponding distribution (in white); mouth corners position for an open mouth.

Figure 8 :

 8 Figure 8: Inner and outer borders detection from the segmented lips and the corresponding mouth corner positions.

Figure 9 :

 9 Figure 9: Convergence snakes results on two sequences (Top: Nico; Bottom: Benny).

Figure 10 :

 10 Figure 10: First results on beard faces sequences: the outer contour.


Mustapha Kchikech
email: kchikech@u-bourgogne.fr

Olivier Togni
email: olivier.togni@u-bourgogne.fr

Approximation Algorithms for Multicoloring Planar Graphs and Powers of Square and Triangular Meshes

Keywords: Graph theory, Coloring, Multicoloring, Planar graph, Power graph, Product graph, Approximation algorithm, Greedy algorithm

A multicoloring of a weighted graph G is an assignment of sets of colors to the vertices of G so that two adjacent vertices receive two disjoint sets of colors. A multicoloring problem on G is to find a multicoloring of G. In particular, we are interested in a minimum multicoloring that uses the least total number of colors.

The main focus of this work is to obtain upper bounds on the weighted chromatic number of some classes of graphs in terms of the weighted clique number. We first propose an 11 6 -approximation algorithm for multicoloring any weighted planar graph. We then study the multicoloring problem on powers of square and triangular meshes. Among other results, we show that the infinite triangular mesh is an induced subgraph of the fourth power of the infinite square mesh and we present 2-approximation algorithms for multicoloring a power square mesh and the second power of a triangular mesh, 3-approximation algorithms for multicoloring powers of semi-toroidal meshes and of triangular meshes and 4-approximation algorithm for multicoloring the power of a toroidal mesh. We also give similar algorithms for the Cartesian product of powers of paths and of cycles.

Introduction

In this paper by graph we mean a simple graph. We denote by G = (V, E) a graph G with vertex set V and edge set E. The length of a path between two vertices is the number of edges on that path. The distance between two vertices is the length of a shortest path between them. We define the Cartesian product G 1 ✷G 2 of two graphs G 1 = (V 1 , E 1 ) and G 2 = (V 2 , E 2 ) to be the graph G = (V, E) with V = V 1 × V 2 and ((x, y), (x ′ , y ′ )) ∈ E whenever x = x ′ and (y, y ′ ) ∈ E 2 or (x, x ′ ) ∈ E 1 and y = y ′ . Given a positive integer p, the p th power G p of a graph G is a graph with the same set of vertices as G and an edge between two vertices if and only if there is a path of length at most p between them in G. A vertex coloring of a graph G is an assignment of colors to the vertices of the graph such that no two adjacent vertices receive the same colors. The minimum number of colors for which a coloring of G exists is called the chromatic number and is denoted by χ(G). Let G = (V, E) be a graph. A weighted graph of G is a pair G ω = (G, ω) where ω is a weight function that assigns non-negative integer ω G (v) (or simply ω(v)) to each vertex v of G, ω G (v) is called the weight of v in G. A vertex multicoloring of the weighted graph G ω consists of a set of colors C and a function Ψ that assigns to each v ∈ V a subset of colors Ψ(v) ⊂ C such that:

i) ∀ v ∈ V , |Ψ(v)| = ω G (v), i-e. the vertex v gets ω G (v) distinct colors.
ii) If (u, v) ∈ E then Ψ(u) ∩ Ψ(v) = ∅, i-e. two adjacent vertices get disjoint sets of colors.

The weighted chromatic number (or multichromatic number), denoted χ ω (G), of G is the minimum number of colors needed to multicolor all vertices of G ω so that conditions i and ii above are satisfied.

A subgraph K of G is called a clique if every pair of vertices in K is connected by an edge. The weight of any clique in G ω is defined as the sum of the weights of the vertices forming that clique. The weighted clique number of G, denoted W (G) (for short, we will use W ), is defined to be the maximum over the weights of all cliques in G ω . Clearly, χ ω (G) ≥ W .

Note that a multicoloring of a weighted graph G = (V, E) is the same as the usual vertex coloring of a graph G ′ that is obtained from G by replacing each vertex u ∈ V by a clique K(u) of size ω(u), and connect two vertices in distinct cliques K(u) and K(v) if and only if u and v are adjacent in G. Nevertheless, depending on the practical problem to be modeled, it is often more convenient to consider it as a multicoloring of a weighted graph since the structure of the graph remains the same even if the weights change. For instance, in the literature, almost all studies about frequency assignment in radio networks adopt this approach.

The multicoloring problem (also known as weighted coloring [START_REF] Mcdiarmid | Channel Assignment and Weighted Coloring[END_REF] or ω-coloring [START_REF] Narayanan | Static Frequency Assignment In Cellular Networks[END_REF]) has been studied in several contexts. Halldórsson and Kortsarz in [START_REF] Halldórsson | Tools for Multicoloring with Application to Planar Graphs and Partial k-Trees[END_REF] present this problem as a model for scheduling of dependent jobs on multiple machines. In [START_REF] Kchikech | Multicoloring and Wavelength Assignment in WDM all-optical Networks[END_REF][START_REF] Nomikos | Routing and path multicoloring[END_REF], a technic based on graph multicoloring has been used to study the routing and wavelength assignment problem on WDM all-optical networks. Also the frequency (channel) assignment problem in a cellular network can be modeled as a multicoloring problem of graph. In the last two decades, this problem was largely studied. Thus, there is a vast literature on algorithms for the multicoloring problem on graphs (especially hexagon graphs), but generally there are no proven bounds on the approximation ratio of the proposed algorithms in terms of the number of colors used in relation to the weighted chromatic number.

The multicoloring problem is NP-hard in general. Hence, it would be interesting to find algorithms that approximate the weighted chromatic number. Note that an algorithm is an f -approximation algorithm for the vertex multicoloring problem if the algorithm runs in polynomial time and it always produces a solution that is within a factor of f of the optimal solution.

Previous results: Exact and approximate algorithms for the multicoloring problem have been given for various classes of graphs, as indicated in Table 1. Note that the multicoloring problem on a triangular mesh has been extensively studied and proved to be NP-hard by McDiarmid and Reed [START_REF] Mcdiarmid | Channel Assignment and Weighted Coloring[END_REF]. This problem corresponds to the frequency (radio channel) assignment in cellular networks. Moreover, if the triangular mesh is of power p ≥ 2 then the problem models channel allocation with interference constraints [START_REF] Feder | Online channel allocation in FDMA networks with reuse constraints[END_REF] also called reuse distance r with r = p -1. Some authors independently gave approximation algorithms for this problem. In case where r = 2, a 4 3 -approximation algorithm has been described both in [START_REF] Mcdiarmid | Channel Assignment and Weighted Coloring[END_REF][START_REF] Narayanan | Static Frequency Assignment In Cellular Networks[END_REF]. For r = 3, [START_REF] Feder | Online channel allocation in FDMA networks with reuse constraints[END_REF] gives a simple algorithm that has a guaranteed approximation ratio of 7 3 . For r ≥ 4, the best known upper bound on the number of colors needed is 4W [START_REF] Jordan | Worst-case preference of cellular channel assignment policies[END_REF], where W is the weighted clique number.

In contrast, the best known lower bound on the number of colors needed is 9 8 W if r = 2 [START_REF] Narayanan | Static Frequency Assignment In Cellular Networks[END_REF] and is 5 4 W if r ≥ 3 [START_REF] Narayanan | Worst-case analysis of a dynamic channel assignment strategy[END_REF].

Our results: In Section 2, we give an algorithm for multicoloring any planar graph with at most 11 6 W colors. In Section 3, we present two upper bounds on the weighted chromatic number of a graph, one in relation with the chromatic number and the other in relation with the minimum number of cliques necessary to cover the neighborhood of any vertex in the graph. Section 4 is devoted to the study of the multicoloring problem on powers of square and triangular meshes. In particular, we show that the infinite triangular mesh is an induced subgraph of the fourth power of the infinite square mesh. We present 2approximation algorithms for multicoloring a power square mesh and the second power of a triangular mesh, 3-approximation algorithms for multicoloring powers of semi-toroidal meshes and of triangular meshes and a 4-approximation algorithm for multicoloring powers of toroidal meshes. In Section 5, we present similar results for graphs of the form G p ✷H q , where G and H are either a path or a cycle.

The best known upper bounds and our results for various classes of graphs are summarized in Table 1. 2 ( * ) (G✷H) p ; G p ✷H q , G is a path; H is a cycle 3 ( * ) p th power toroidal mesh 4 ( * )

Graphs

Tab. 1: Best approximation ratios for multicoloring algorithms of some classes of graphs. ( * ) our results in this paper.

Multicoloring planar graphs

In this section we consider the multicoloring problem on planar graphs. Before presenting our result, we remind a result of Narayanan and Shende [START_REF] Narayanan | Static Frequency Assignment In Cellular Networks[END_REF] showing that there exists an efficient algorithm to optimally multicolor any outerplanar graph. A graph is planar if it can be drawn in a plane without edge crossings.

A graph is said to be outerplanar if it is a planar graph so that all vertices may lie on the outer face.

Theorem 1 ([12]

) Let G be an arbitrary outerplanar graph, then its associated weighted graph G ω can be multicolored optimally using χ ω (G) colors in linear time.

Theorem 2 Let G be a planar graph, then

χ ω (G) ≤ 11 6 W (G) .
Proof: We may assume G = (V G , E G ) to be connected, since disconnected components of G can be multicolored independently. As G is a planar graph of order n, then using the O(n 2 ) algorithm described by Robertson et al. in [START_REF] Robertson | The four-colour theorem[END_REF], we color G with 4 colors from {1, 2, 3, 4}. We call these colors base colors. We denote by s i any vertex s ∈ V G which has color i ∈ {1, 2, 3, 4}, and by {c 1 i , c 2 i , • • • , c α i } a set of α (nonnegative integer) distinct hues associated with base color i so that if i = j then for every integers p, q ≥ 1, we have c p i = c q j .

Case 1: W (G) ≡ 0 (mod3).

We fix ℓ = 1 3 W (G), and we let

C = 4 i=1 {c 1 i , c 2 i , • • • , c ℓ i } denote a set of available colors.
Consider the multicoloring function f of G ω defined as follows:

f : V -→ P(C) s i -→ f (s i ) = {c 1 i , c 2 i , • • • , c α i }, with α = min(ω G (s i ), ℓ)
where each vertex

s i with weight ω G (s i ) is assigned the hues of {c 1 i , c 2 i , • • • , c α i }. We call a vertex v heavy if ω(v) > ℓ or light if ω(v) ≤ ℓ.
Hence, only the heavy vertices remain to be completely colored and their weights may be decreased by ℓ. All light vertices are colored completely and are deleted from G.

Let H = (V H , E H ) denote the remaining graph obtained after this process. Thus H is such that

• u ∈ V H ⇔ ω G (u) > ℓ, • ω H (u) = ω G (u) -ℓ, • (u, v) ∈ E H ⇒ (u, v) ∈ E G .
It is easily seen that H has no clique of size 3. In fact, if (x, y, z) is a triangle in H, then these vertices must have been heavy in G. Hence, there exist positive integers ε 1 , ε 2 and

ε 3 such that ω G (x) = ℓ + ε 1 , ω G (y) = ℓ + ε 2 and ω G (z) = ℓ + ε 3 . As ω G (x) + ω G (y) + ω G (z) ≤ W G = 3ℓ, we obtain ε 1 + ε 2 + ε 3 ≤ 0 a contradiction. Consequently, every clique K in H has size at most 2. Furthermore, if (u, v) ∈ K then ω H (u) + ω H (v) = ω G (u) -ℓ + ω G (v) -ℓ ≤ W G -2ℓ = ℓ.
In addition, if u is an isolated vertex in H, then there exists a positive integer ε such that ω G (u) = ℓ + ε, and all neighbors of u must be light vertices in G. Suppose that u has i for its base color and that it has a neighbor in G of base color j = i. Let N j be the set of all its neighbors in G of base color j.

Let v j ∈ N j such that ω G (v j ) = max s∈Nj ω G (s). Further, we have ω G (v j ) = ℓ -ε vj with 0 ≤ ε vj ≤ ℓ since u is isolated in H. Then u can borrow ε vj available colors from {c 1 j , c 2 j , • • • , c ℓ j }, which are unused by all vertices of N j . As ω G (u) + ω G (v j ) ≤ W G = 3ℓ, we get ε -ε vj ≤ ℓ. For this reason, we consider that each isolated vertex u ∈ H has ω H (u) = ω G (u) -(ℓ + ε vj ) = ε -ε vj .
Thus, W H ≤ ℓ and we can therefore distinguish two cases: (a) If H is a bipartite graph, then H can be multicolored optimally with exactly W H colors (see [START_REF] Janssen | Fixed preference frequency allocation for cellular telephone systems[END_REF][START_REF] Mcdiarmid | Channel Assignment and Weighted Coloring[END_REF]). In this case, to avoid color conflicts, we use a new set C ′ of W H distinct colors. Thus, multicoloring all vertices of G ω requires at most

|C| + |C ′ | ≤ 4ℓ + ℓ = 5 3 W G colors.
(b) If H is not a bipartite graph, then as H is a planar graph without 3-cycles (triangle-free), it is well known that H is 3-colorable [START_REF] Grötzsch | Ein dreifarbensatz fur dreikreisfreie netze auf der kuzel[END_REF]. Using the O(n log n) algorithm of Kowalik [START_REF] Kowalik | Fast 3-Coloring Triangle-Free Planar Graphs[END_REF], we color H with 3 colors. Thus, there is an algorithm for multicoloring H that requires at most 3 2 W H colors (See Theorem 4 or [START_REF] Janssen | Fixed preference frequency allocation for cellular telephone systems[END_REF]). Similarly, we use a new set C ′ of 3 2 W H distinct colors. Thus, multicoloring all vertices of G ω requires at most

|C|+|C ′ | ≤ 4ℓ+ 3 2 ℓ = 11 6 W G colors. Case 2: W (G) ≡ 1 (mod3).
In this case, the main idea remains the same as that of Case 1, except that we fix four integers

ℓ 1 , ℓ 2 , ℓ 3 , ℓ 4 such that ℓ 1 = W +2 3 and ℓ 2 = ℓ 3 = ℓ 4 = W -1 3 .
Note that ℓ i represents the number of distinct hues associated with base color i. Then, each vertex s i of base color i is assigned the hues of

{c 1 i , c 2 i , • • • , c α i }, with α = min(ω(s i ), ℓ i ). A vertex v i is called heavy if ω(v i ) > ℓ i and is called light if ω(v i ) ≤ ℓ i .
All the light vertices thus get completely colored and are deleted from G. The weight of every remaining heavy vertex u i is decreased by ℓ i , resulting in the graph H.

Observe that every triangle in G must contain at least one light vertex. Otherwise, there exist three heavy vertices x, y, z of G such that (x, y, z) is a triangle in H. Hence, there exist positive integers

ε 1 , ε 2 , ε 3 ≥ 1 such that ω G (x i ) = ℓ i + ε 1 , ω G (y j ) = ℓ j + ε 2 and ω G (z k ) = ℓ k + ε 3 . As ω G (x i ) + ω G (y j ) + ω G (z k ) ≤ W G , we obtain ε 1 + ε 2 + ε 3 ≤ 0 or ε 1 + ε 2 + ε 3 ≤ 1 (because ℓ i + ℓ j + ℓ k = W G -1 if i, j, k ∈ {2, 3, 4} or ℓ i + ℓ j + ℓ k = W G otherwise) a contradiction.
Consequently, H has no cliques of size greater than 2.

Let K be a clique of H. If |K| = 2 then K contains two vertices u i and v j such that ω H (u i ) + ω H (v j ) = ω G (u i ) -ℓ i + ω G (v j ) -ℓ j ≤ W G -(ℓ i + ℓ j ) ≤ ℓ 1 . If |K| = 1 then K contains an isolated vertex u i such that ω G (u i ) = ℓ i + ε, with ε ≥ 1,
and all neighbors of u i must be light vertices in G. Assume u i has a neighbor in G of base color j = i and let N j (u i ) be the set of all neighbors of u i in G of base color j.

Let v j ∈ N j (u i ) such that ω G (v j ) = max s∈Nj (ui) ω G (s).
Further, we have ω G (v j ) = ℓ jε vj with 0 ≤ ε vj ≤ ℓ j since u i is isolated in H. Then u i can borrow ε vj available colors from {c 1 j , c 2 j , • • • , c ℓj j }, which are unused by all vertices of N j (u i ).

As ω G (u i ) + ω G (v j ) ≤ W G , we get ε -ε vj ≤ W G -(ℓ i + ℓ j ) ≤ ℓ 1 (because ℓ i + ℓ j = 2W G -2 3 if i, j ∈ {2, 3, 4} or ℓ i + ℓ j = 2W G +1 3 otherwise). For this reason, we consider that each isolated vertex u i ∈ H has ω H (u i ) = ω G (u i ) -(ℓ i + ε vj ) = ε -ε vj ≤ ℓ 1 .
Consequently, for any clique K of H we have ω H (K) ≤ ℓ 1 . Thus, W H ≤ ℓ 1 and so H can be multicolored using a new set C ′ of colors, with

|C ′ | = ℓ 1 if H is bipartite or |C ′ | = 3
2 ℓ 1 otherwise. Then the total number of colors used in this case is

|C|+|C ′ | ≤ ℓ 1 +ℓ 2 +ℓ 3 +ℓ 4 + 3 2 ℓ 1 = 11 6 W G + 2 3 . Case 3: W (G) ≡ 2 (mod3).
Let ℓ 1 , ℓ 2 , ℓ 3 , ℓ 4 be four non-negative integers such that ℓ 1 = W -2

3
and ℓ 2 = ℓ 3 = ℓ 4 = W +1 3 , and we can use the same method we used in Case 2.

✷ 3 Relationships between W (G) and χ ω (G)

The problem of precisely determining the weighted chromatic number is an intractable one in general. Nonetheless, it is possible to determine the exact value of the weighted chromatic number on specific families of graphs, such as bipartite graphs [START_REF] Mcdiarmid | Channel Assignment and Weighted Coloring[END_REF] and outerplanar graphs [START_REF] Narayanan | Static Frequency Assignment In Cellular Networks[END_REF]. In this section, we provide two upper bounds for χ ω (G), one in relation with the chromatic number and the other in relation with the minimum number of cliques necessary to cover the neighborhood of any vertex in the graph. The second result plays a role in the rest of our work where it helps to determine upper bounds for weighted chromatic numbers of various classes of graphs of the form (G✷H) p or G p ✷H q .

Given a k-colorable graph G (then G have chromatic number χ(G) ≤ k), Janssen, Kilakos and Marcotte [START_REF] Janssen | Fixed preference frequency allocation for cellular telephone systems[END_REF] showed that there exists a k 2 -approximation algorithm for multicoloring all vertices of G.

Theorem 3 ( [5]) If G is a k-colorable graph, then χ ω (G) ≤ k 2 W (G).
However, finding a maximal clique has been proven to be an N P -complete problem. Thus, determining W (G) becomes also a difficult problem. So, to calculate the upper bound on χ ω (G) by using Theorem 3 is complicated. In this section, we present an improvement of the above result using an approximation algorithm that can be implemented in a distributed manner.

In the following, we denote by W t (G) = max{W (K) | K clique of G s.t. |K| ≤ t} the maximum over the weights of all cliques of size at most t in G ω .

Theorem 4 If G is a k-colorable graph without isolated vertices, then

χ ω (G) ≤ k 2 W 2 (G).
Proof: The proof is based on the same general principle as that of Theorem 2. We may assume G to be connected, since disconnected components of G can be multicolored independently. As G is a kcolorable graph, then we color G with at most k colors from {1, 2, ..., k}. We call these colors base colors. We denote by v i any vertex v of G which has base color i, and by {c 1 i , c 2 i , • • • , c α i } a set of α (nonnegative integer) distinct hues associated with base color i so that if i = j then for every integers p, q ≥ 1, we have c p i = c q i . First, we assume W 2 (G) ≡ 0 (mod2) and we fix ℓ

= 1 2 W 2 (G) (if W 2 (G) ≡ 1 (mod2)
, we can use similar arguments as those presented in the proof of Theorem 2). Every vertex v i is assigned the first ω(v i ) hues of the set

{c 1 i , c 2 i , • • • , c ℓ i }.
Hence, only the heavy vertices (of weight at least ℓ + 1) remain to be completely colored after this step and the other vertices are completely colored and are deleted from G. Let H denote the remaining graph obtained after this process where the new weight of each vertex v ∈ H is ω H (v) = ω(v)ℓ. We observe that H contains only isolated vertices, because if we suppose that there exist two heavy adjacent vertices v and u in H then, we get 2ℓ < ω(v) + ω(u) ≤ W 2 (G) = 2ℓ, a contradiction. Assume v j has a neighbor in G of base color i = j and let N i (v j ) be the set of neighbors vertices of v j having i as base color in G.

Let v i ∈ N i (v j ) such that ω(v i ) = max v∈Ni(vj ) ω(v). Further, we have ω(v j ) = ℓ + ǫ j and ω(v i ) = ℓ -ǫ i where ǫ j > 0 and 0 ≤ ǫ i ≤ ℓ. As ω(v j ) + ω(v i ) ≤ W 2 (G) = 2ℓ,
we get ǫ j ≤ ǫ i . Then v j can borrow from ǫ i colors available in the set {c 1 i , c 2 i , • • • , c ℓ i } for coloring the remaining weight on v j . Consequently, for multicoloring all vertices of G, we use at most kℓ

= k 2 W 2 (G) colors. Thus, χ ω (G) ≤ k 2 W 2 (G). ✷
Remark that, for a graph G on m edges, W 2 (G) can be computed in time O(m) and W 2 (G) ≤ W (G), thus our result is an improvement of Theorem 3. Now, we give a upper bound on the ratio χω(G) W (G) in relation with the minimum number of cliques necessary to cover the neighborhood of any vertex.

For a vertex u of a graph

G = (V, E), let N G [u] be the closed neighborhood of u: N G [u] = {v ∈ V |u = v or uv ∈ E}.
Definition 1 A connected graph G has Property P q if and only if for any vertex u, N G [u] can be covered with q distinct cliques each containing u.

Lemma 1 Let G = (V, E) be a simple graph and ω be a weighting.

1. If q ∈ N is such that G has Property P q , then the vertices of G can be multicolored in a greedy fashion using at most q(W -1)+ 1 colors in any given order. In particular χ ω (G) ≤ q(W -1)+ 1.

2. If further there is a q ′ ≤ q such that for a specific ordering

V = {u 1 , u 2 , • • • , } (not necessarily finite) of the vertices of G, we have that N G [u i ] ∩ {u 1 , • • • , u i-1
} can be covered by at most q ′ cliques, then the vertices can be multicolored in a greedy fashion using at most q ′ W colors in this order. In particular χ ω (G) ≤ q ′ (W -1) + 1.

Proof:

1. Assume G is k-colorable and we associate with each base color ℓ distinct hues, where ℓ is some fraction of W to be determined later. Thus, we obtain kℓ available colors. Now, we assign to each vertex u i of G of color i α = min(ω G (u i ), ℓ) colors in the set {c 1 i , c 2 i , • • • , c ℓ i }, associated to the base color i of u. Then, we note that the heavy vertices (of weight at least ℓ + 1) are not completely colored. In order to complete the coloring of the heavy vertices, we proceed in the following manner.

Consider a heavy vertex u not yet completely colored. As G satisfies Property P q , N G [u] can be covered with q cliques, each of them containing u. Note that each clique has weight at most W , since W is the weighted clique number. Thus, the total weight of

N G [u] is at most ω(u) + q(W - ω(u)) = qW -(q -1)ω(u) ≤ qW -(q -1)(ℓ + 1) because ω(u) ≥ ℓ + 1.
Consequently, kℓ colors will be sufficient to color all vertices of N G [u] if qW -(q-1)(ℓ+1)ℓ ≤ kℓ. This gives q(W -1) ≤ (q -1 + k)ℓ -1. So, we can take ℓ = q(W -1)+1 q-1+k . Thus, we are able to color all vertices of G using at most kℓ = k(q(W -1)+1) q-1+k ≤ q(W -1) + 1 colors.

We consider the vertices of

G in a specific order V (G) = {u 1 , u 2 , • • •} such that N G [u i ] ∩ {u 1 , • • • , u i-1
} can be covered by at most q ′ cliques. The idea of the proof is to multicolor in a greedy fashion each vertex u i of G using the same method we used in the first assertion. approximation algorithm for multicoloring a finite square mesh with 4 3 W + 4 colors. In the following three theorems, we present algorithms for multicoloring powers of square meshes, semi-toroidal meshes and toroidal meshes, with approximation ratios of 2, 3 and 4 respectively. Their proofs are direct consequences of Lemma 1 and Lemma 2.

Theorem 5 For any p ≥ 2, there exists a polynomial time greedy algorithm which multicolors all vertices every finite induced portion of a weighted p th power of the infinite square mesh M p using at most (2W -1) colors.

Theorem 6 Consider (G✷H) p , where one of the graphs G or H is a cycle. There exists a polynomial time greedy algorithm which multicolor all vertices of every finite induced portion of (G✷H) p using at most (3W -2) colors.

Theorem 7 Consider the p th power toroidal mesh T M p = (G✷H) p where both G and H are cycles. There exists a polynomial time algorithm which multicolors all vertices of T M p with using at most (4W -3) colors.

Powers of triangular meshes

Now, we consider the multicoloring problem on the p th power of a triangular mesh. We define the infinite triangular mesh H as a mesh formed by tiling the plane regularly with equilateral triangles.

In the following, we present an approximation algorithm that multicolor all vertices of any finite induced portion of H p using at most (2W -1) colors if p = 2 and at most (3W -2) colors if p ≥ 3. The method used is based on the multicoloring of the p th power of a square mesh.

Theorem 8

The infinite triangular mesh H is an induced subgraph of (Z✷Z) 4 , the fourth power of the infinite square mesh.

Proof: The infinite triangular mesh can be presented by H = {x(1, 3) + y(3, 1) : x, y ∈ Z} = {(x + 3y, 3x + y) : x, y ∈ Z} ⊆ Z 2 , where the edges are between all pairs of points of distance four apart (See Fig. 4). ✷ 00000000000000000000000 00000000000000000000000 00000000000000000000000 00000000000000000000000 00000000000000000000000 00000000000000000000000 00000000000000000000000 00000000000000000000000 00000000000000000000000 00000000000000000000000 00000000000000000000000 00000000000000000000000 00000000000000000000000 00000000000000000000000 As a consequence of Theorem 8, for every finite induced portion Γ of H, there are m and n such that Γ embeds into (P m ✷P n ) 4 .

Theorem 9 There exists a polynomial time greedy algorithm which multicolor all vertices of every finite induced portion of the 2 th power of a triangular mesh H 2 using at most (2W -1) colors.

Proof: If we let H 2 be the subgraph of M 8 defined previously (see proof of Theorem 8), then we can use the same multicoloring process used in the proof of Theorem 5. In fact, keeping the same order on the vertices, when we multicolor vertex u k of H 2 , the previously listed neighbors of u k can be covered by two cliques. Thus the same algorithm used in proof of Theorem 5 gives the result. ✷ Theorem 10 For any p ≥ 3, there exists a polynomial time greedy algorithm which multicolor all vertices of every finite induced portion of the p th power of a triangular mesh H p using at most (3W -2) colors.

Proof: The proof is similar to that of Theorem 5. However, when we multicolor the vertex u k of H p , one can see that the previously listed neighbors of u k can be covered by three cliques. ✷

Multicoloring the product of powers of paths and of cycles

Consider G p ✷H q where both G and H are either a simple path (finite or infinite) or a simple cycle. The multicoloring problem on this type of graphs is of interest since it can model the problem of coloring a set of line-column paths in a mesh [START_REF] Kchikech | Paths coloring algorithms in mesh networks[END_REF][START_REF] Kchikech | Multicoloring and Wavelength Assignment in WDM all-optical Networks[END_REF]. Note that G p ✷H q is a subgraph of (G✷H) r , with r = max(p, q). But, as it is not an induced subgraph, the results for powers of square meshes cannot be extended directly to graphs of the form G p ✷H q . Nevertheless, the analogous of Lemma 2 for graphs G p ✷H q is easy to state. Lemma 3 Consider G p ✷H q where both G and H are either a simple path (finite or infinite) or a simple cycle.

1. Any closed neighborhood of a vertex can be covered by four cliques of G p ✷H q .

2. If exactly one of the graphs G or H is a cycle, then there is an ordering of the vertices, such that the previously listed neighbors of each vertex can be covered by three cliques of G p ✷H q .

3. If neither of the graphs G nor H are cycles, then there is an ordering of the vertices, such that the previously listed neighbors of each vertex can be covered by two cliques of G p ✷H q .

Notice that, in the case where both G and H are a simple (finite or infinite) path, it is known that χ(G p ✷H q ) = max(χ(G p ), χ(H p )) = max(p + 1, q + 1). Thus, for the vertex multicoloring problem on G p ✷H q , if we apply Theorem 3, we obtain χ ω (G p ✷H q ) ≤ max(p+1,q+1) 2 W . Consequently, we have χ ω (G p ✷H q ) ≤ 3 2 W for the case where max(p, q) ≤ 2. General upper bounds on χ ω (G p ✷H q ) are given in the following theorem whose proof is a direct consequence of Lemma 1 and Lemma 3.

Theorem 11 Let Γ be any finite induced portion of G p ✷H q , where both G and H are either a simple path (finite or infinite) or a simple cycle. There exists a polynomial time greedy algorithm which multicolor all vertices of Γ, using a total of at most

Fig. 3 :

 3 Fig.3: The graph (P8✷P10)3 with the previously listed neighbors of ui in deglex order being covered by 2 cliques (only edges of (P8✷P10) 1 are represented).

Fig. 4 :

 4 Fig.4: The infinite triangular mesh induced in the 4 th power of the infinite square mesh.

  1. (2W -1) colors if neither of the graphs G nor H are cycles, 2. (3W -2) colors if exactly one of the graphs G or H is a cycle, 3. (4W -3) colors if both G and H are cycles.

1365-8050 c 2006 Discrete Mathematics and Theoretical Computer Science (DMTCS), Nancy, France

Acknowledgements

The authors would like to thank the anonymous referees for many helpful comments which greatly improved the presentation of the paper.

Multicoloring powers of square and triangular meshes

In this section we present greedy algorithms for multicoloring square and triangular meshes and we compute their approximation ratios.

Powers of square meshes

We begin with the following easy to observe lemma.

Lemma 2 Consider (G✷H) p , where both G and H are either a simple path (finite or infinite) or a simple cycle.

1. Any closed neighborhood of a vertex can be covered by four cliques of (G✷H) p .

2. If exactly one of the graphs G or H is a cycle, then there is an ordering of the vertices, such that the previously listed neighbors of each vertex can be covered by three cliques of (G✷H) p .

3. If neither of the graphs G nor H are cycles, then there is an ordering of the vertices, such that the previously listed neighbors of each vertex can be covered by two cliques of (G✷H) p .

Proof: Without loss of generality, we denote by (i, j), with (i, j) ∈ N 2 , the vertices of (G✷H) p .

1. Let s = (i, j) be a vertex of (G✷H) p , and let N (G✷H) p [s] be the closed neighborhood of s.

Consider the following four subsets of V ((G✷H) p ):

It is easy to see that each subset K l ij with l ∈ {1, 2, 3, 4} is a clique in V ((G✷H) p ) and that N (G✷H) p [s] is covered by these four cliques. An illustration is given in Fig. 1 for the case p = 3 (only the edges of (G✷H) 1 are drawn for clarity).

2. If exactly one of the graphs G or H is a cycle, then the vertices of (G✷H) p are ordered using the degree lexicographical order (deglex for short) [START_REF] Adams | An introduction to Gröbner bases[END_REF] of the points of

It is then easy to observe that for each vertex

} can be covered by three cliques. As above, an illustration is given in Fig. 2 for the case p = 3 (only the edges of (G✷H) 1 ) are drawn for clarity).

3. If neither of the graphs G nor H are cycles, then, as previous, listing the vertices of (G✷H) p in deglex order yelds the desired condition that for each vertex

can be covered by two cliques (See Fig. 3).

Fig. 1: The graph (P8✷P10) 3 with the neighbors of (i, j) being covered by

✷

Let M be the infinite square mesh with vertex set Z✷Z, where two vertices (x 1 , y 1 ) and (x 2 , y 2 ) of M are adjacent if and only if

Consider the infinite weighted p th power mesh M p with weighted clique number W . If p = 1 then M 1 is a bipartite graph. McDiarmid and Reed in [START_REF] Mcdiarmid | Channel Assignment and Weighted Coloring[END_REF] proposed an algorithm to optimally multicolor this family of graphs. In addition, when p = 2, the multicoloring of a subgraph of M 2 called lattice graph with diagonals (i) has been studied by Miyamato and Matsui in [START_REF] Miyamoto | Linear Time Approximation Algorithm for Multicoloring Lattice Graphs with Diagonals[END_REF]. They proposed an O(|M |) time (i) Graph obtained from the standard square mesh by adding the two diagonals edges on each face.


Pierrick Coupé

Olivier Périot

Jose Vicente Manjon

Bassem Hiba

Michèle Allard

Olivier Periot

José V Manjón

Impact of DWI denoising on Track-Density Imaging

  L'archive ouverte pluridisciplinaire

Introduction

Track-Density Imaging (TDI) [START_REF] Calamante | Track density imaging (TDI): Validation of super resolution property[END_REF] has been proposed as super-resolution method for diffusion-weighted images (DWI). Based on fiber-tracking (FT), TDI enables to produce high-resolution white matter images. However, the sensitivity to noise of FT negatively impacts TDI accuracy and reproducibility [START_REF] Dhollander | Track-density imaging + noise: when super-resolution quality does not yield accuracy[END_REF]. In this study, we analyzed whether the use of DWI denoising is able to improve TDI robustness and reproducibility.

Methods

Data: 10 DWI acquisitions of a human brain were acquired using a 3T MR scanner (Achieva, Philips Medical Systems). A standard spin-echo EPI pulse sequence with sensitivity encoding was used (b-value: 700 s/mm2; 21 diffusion directions) with the following parameters: TE/TR ~ 60/8000 ms; FOV = 215x215x85 mm 3 ; matrix size = 172x172 with 68 slices and an isotropic spatial resolution of 1.2x1.2x1.2 mm 3 .

Image Processing:

To compensate for inter-acquisition motion and intra-acquisition eddy current distortion, all the DWI (i.e., 220 images) were linearly registered to the b0 image of the first acquisition. As shown in [START_REF] Coupé | Collaborative patchbased super-resolution for diffusion-weighted images[END_REF], trilinear interpolation of FSL eddy current correction tends to overblur the corrected images. Therefore, the FSL script was modified to use nearest neighbor interpolation to better preserve structures and noise patterns. Finally, the b-matrixes of the 10 acquisitions were adjusted to account for head rotation.

Noisy and Denoised testing datasets:

To investigate the impact of denoising on TDI, two testing datasets were used. First, the noisy dataset was composed of the 10 motion-corrected original acquisitions. Second, the denoised dataset was based on the same 10 acquisitions but denoised using the LPCA denoising filter [START_REF] Manjon | Diffusion Weighted Image Denoising using overcomplete Local PCA[END_REF].

Gold-Standard (GS):

The 10 acquisitions of the noisy dataset were averaged in order to estimate the GS. This resulted in a high SNR GS used as reference during experiments. TDI: TDI [START_REF] Calamante | Track density imaging (TDI): Validation of super resolution property[END_REF] was estimated from the GS, the noisy and the denoised datasets. The diffusion tensors and FA maps were obtained by using FSL. Whole brain tractography was obtained with MRtrix by using constrained spherical deconvolution [START_REF] Tournier | Robust determination of the fibre orientation distribution in diffusion MRI: non-negativity constrained super-resolved spherical deconvolution[END_REF] and probabilistic streamlines tractography. The default tracking parameters were used (seed = whole brain, step-size = 0.2 mm, maximum harmonics order = 6, 10 000 000 tracks). Finally, a super-resolution TDI map was created at 0.3 mm isotropic resolution.

Results

To investigate the impact of LPCA denoising on TDI, PSNR was estimated between TDI of both testing datasets and TDI of the GS. This quality metrics was computed over a mask where TDI of the GS was higher than 0. Higher values of PSNR indicate higher similarity compared to GS. The average PSNR was 24.14dB (+/-0.28) for the noisy dataset and 27.30dB (+/-0.14) for the denoised dataset. Figure 1 shows TDI maps calculated from GS, noisy and denoised datasets. The structures are visually more distinct and less noisy in the denoised TDI map. Figure 2 shows the maps of coefficient of variation (CoV)=100*σ/µ and of coefficient of dispersion (CoD)=σ 2 /µ estimated over the 10 acquisitions for both datasets. CoV and CoD maps of the noisy dataset demonstrate a more pronounced speckled appearance with high values in GM and at the level of the ventricles.

Conclusion

The gain of 3dB in term of PSNR demonstrates the higher similarity between TDI of denoised acquisitions and TDI of the GS. Moreover, the lower standard deviation of this metric highlights the better consistency of TDI estimation after denoising. Visual assessment of RGB TDI maps demonstrates a better contrast and conspicuity of boundaries between GM, WM and ventricles on the denoised images. Moreover, in CoV and CoD maps, high values indicate regions where TDI measurements are the less consistent over the 10 acquisitions. By using LPCA filter, it appears that TDI reproducibility is globally improved. Finally, DWI denoising enables to limit false positive TDI values in GM and CSF.


Vincent Levorato
email: vlevorato@cesi.fr

Core decomposition in Directed Networks: Kernelization and Strong Connectivity

Keywords: Kernelization and Strong Connectivity cores, communities detection, directed networks, kernelization, strong connectivity, graph theory, clustering

Introduction

Complex networks appear in many applications, including social networks analysis on the Web, which is a topical research subject. These networks carry non-trivial topological properties that characterize their connectivity, and affect the dynamics of their behaviors. The analysis of complex networks often leads to the analysis of the roles of elements, or groups of elements, composing a network. Communities detection belongs to this research field, and can be very useful to better understand how networks are structured. In this article, we focus on the problem of finding communities in networks, and more specifically finding cores in directed networks. Dealing with methods of community detection for directed networks is a difficult task, and few methods exist compared to methods used in the undirected case. Here are some of the most known works in the literature [START_REF] Fortunato | Community detection in graphs[END_REF][START_REF] Lancichinetti | Community detection algorithms: A comparative analysis[END_REF][START_REF] Labatut | Detection and interpretation of communities in complex networks: Practical methods and application[END_REF] dedicated to directed networks, or which can be adapted to work with directiveness: Clauset et al. method [START_REF] Clauset | Finding community structure in very large networks[END_REF], CFinder based on the Clique Percolation Method (CPM) [START_REF] Palla | Uncovering the overlapping community structure of complex networks in nature and society[END_REF], Louvain method [START_REF] Vincent | Fast unfolding of communities in large networks[END_REF], InfoMap [START_REF] Rosvall | Maps of random walks on complex networks reveal community structure[END_REF], Simulated Annealing for modularity [START_REF] Guimerà | Functional cartography of complex metabolic networks[END_REF], Wu-Huberman method [START_REF] Wu | Finding communities in linear time : A physics approach[END_REF], MarkovCluster algorithm [START_REF] Van Dongen | Graph clustering via a discrete uncoupling process[END_REF], Multistep Greedy algorithm [START_REF] Schuetz | Efficient modularity optimization by multistep greedy algorithm and vertex mover refinement[END_REF] and EM method (Expectation-Maximization) [START_REF] Leicht | Community structure in directed networks[END_REF]. More recently, others methods concerning directed networks have been proposed, with more or less good results [START_REF] Yang | Directed network community detection: A popularity and productivity link model[END_REF][START_REF] Lambiotte | Flow graphs: Interweaving dynamics and structure[END_REF][START_REF] Lai | Finding communities in directed networks by pagerank random walk induced network embedding[END_REF]. Generally, these methods are most of the time designed for undirected networks, and adjusted to work in the directed case: they are not initially dedicated to the directed case. There is also a significant amount of methods using modularity optimization. However, this kind of approach has its limits, and can "miss important substructures of a network" [START_REF] Fortunato | Resolution limit in community detection[END_REF]. Some recent work discuss about reciprocated interaction [START_REF] Cheng | Predicting reciprocity in social networks[END_REF], that two people should communicate in both directions, the first person expecting messages from the second person, and vice versa. Our approach is based on this simple idea: in a directed network, a community should be composed of nodes which can communicate with every nodes in the community, in both directions. Interesting results in our previous exploratory work [START_REF] Levorato | Detection of communities in directed networks based on strongly p-connected components[END_REF] encourage us to continue in this direction. Usually represented by graphs in undirected networks, this kind of representation can be modeled by the connected component concept, and more restrictively by the clique concept. Except that for the directed case, it can be represented by the concept of strongly connected component. Finding these components should be equivalent to find cores to which other elements of the network will be assigned. This work gives the key concepts of this approach, focusing on the cores finding. Our paper is structured as such: the first section gives the definitions of graph theory and formal concepts needed to understand our method. Then, the second part exposes the different steps of our approach, followed by some experimental results on real networks. The paper ends by a conclusion which opens discussion on future work directions.

Graph Theory Notions

Graph definitions

In this article, we consider only directed graphs (also noted digraph). We give here a short reminder of graph theory notions. Formally, a digraph G = (V, A) is the pair composed of [START_REF] Berge | Graphes et Hypergraphes[END_REF]:

-a set V = {x 1 , x 2 , ..., x n } named vertices or nodes. -a family A = (a 1 , a 2 , ..., a n ) of elements of the Cartesian product V × V = {(x, y)/x ∈ V, y ∈ V } named arcs.
The amount of vertices is noted n (also noted |V (G)|) and the amount of arcs is noted m (also noted |A(G)|). A path P is composed of k arcs such as P = (a 1 , a 2 , ..., a i , ..., a k ) where for every arc a i the terminal end coincides with the initial end of a i+1 . Several equivalent notations can be used:

P = ((x 1 , x 2 ), (x 2 , x 3 ), ...) = [x 1 , x 2 , ..., x k , x k+1 ] = P [x 1 , x k+1 ].
A chain is, like a path, an alternating sequence of vertices and edges, where an edge is an arc without orientation. A circuit is a path such that the first node of the path corresponds to the last. It can be viewed as an oriented cycle.

Connected Components

Here are the different types of connected components we could have in a directed graph [START_REF] Harary | Graph Theory[END_REF]:

a weakly connected component W CC of a digraph is a subgraph where: ∀x, y ∈ W CC, there is a chain between x and y. an unilaterally connected component U CC of a digraph is a subgraph where: ∀x, y ∈ U CC, there is a path between x and y OR there is a path between y and x. a strongly connected component SCC of a digraph is a subgraph where: ∀x, y ∈ SCC, there is a path between x and y AND there is a path between y and x.

To discover cores in a network, we use a special case of strongly connected component named strongly p-connected component by [START_REF] Wasserman | Social Network Analysis: Methods and Applications[END_REF] which is related to l-edgeconnectivity [START_REF] Diestel | Graph Theory[END_REF] (we use p-connected notation instead of n-connected to avoid confusion):

a strongly p-connected component p-SCC of a digraph is a subgraph where: ∀x, y ∈ p-SCC, there is a path of length p or less between x and y, and there is a path of length p or less between y and x, with p ≥ 2.

3 Core detection

Related work

Finding cores in order to find communities is a method that can be related to pattern identification [START_REF] Labatut | Detection and interpretation of communities in complex networks: Practical methods and application[END_REF]. This consists in finding maximal subsets which implies separation between them. Clique finding is one of these methods, but is also very restrictive, because each node must have a direct connection to other nodes. This approach has been relaxed by the n-clique definition where each node is connected to others by at least one path which length is at most n, but that can be outside of the n-clique. The n-clan concept, by adding a constraint on the diameter of a n-clique that should not be greater than n, fixes the connectedness issue of the n-clique [START_REF] Mokken | Cliques, clubs and clans[END_REF]. Some concepts in the directed case exist such as the f -groups (maximal subsets of weakly and strongly transitive triads, like 3-clique or triangle) [START_REF] Hanneman | Introduction to social network methods[END_REF], and directed k-clique [START_REF] Palla | Directed network modules[END_REF]. Our approach is related to these works, but we don't put strong constraint on the size (triads), and we don't want to avoid circuits (directed k-clique), as it is specifically the configuration we are looking for: strongly p-connected components.

Searching for cores using p-SCC concept

Our community definition refers to a group containing elements that can communicate with all other elements of the group. In digraphs, this idea is formalized by the concept of strongly connected components (SCC). Searching for SCCs in a digraph is equivalent to search for circuits, and Tarjan based his algorithm on this idea [START_REF] Tarjan | Depth-first search and linear graph algorithms[END_REF], which returns maximal strongly connected subgraphs. To our knowledge, no work has considered the SCCs in the case of researching communities. By simply applying Tarjan's algorithm on directed graph generated through LFR benchmark [START_REF] Lancichinetti | Benchmarks for testing community detection algorithms on directed and weighted graphs with overlapping communities[END_REF], some communities can be found, but SCCs are often oversized. To refine the process, our approach proposes to find p-SCCs (fig. 1). The problem is that in a digraph, the number of circuits may be exponential in the number of vertices [START_REF] Tarjan | Enumeration of the Elementary Circuits of a Directed Graph[END_REF]. Therefore, processing all circuits of a graph is not relevant, especially if the graph has a significant number of nodes like in large real-world networks.

(a) (b) To be able to process large graph, our method finds p-SCCs by starting from a given node s. Starting from this node, finding p-SCC also means searching for circuits, but circuits with a given size. Trivially, the length of the path p is bounded by the length of the circuits found into the p-SCC :

p ≤ 2 × (c -1)
where c is the size of circuits we are searching for. s Fig. 2. Searching p-SCC is similar to search circuits from a starting node s. In this case, searching for 4-SCC means searching for circuits of length at most 3. The highlighted path length is 4, the maximal path length which can be found.

For instance, searching for circuits of length 3 starting from a node means searching for at most 4-SCCs (fig. 2). We propose an algorithm which returns a p-SCC starting from a given node (alg. 1). As the algorithm is searching for circuits, and considering that p parameter sets the circuit length, we can only find p-SCCs with p being an even number. It is written in a non-recursive way, but time complexity should be approximatively the same as Tarjan's algorithm, which is O(n + m), with two differences: we don't always need to pass through every arc (depends on path length), but we should pass through nodes several times. Algorithm 1: Algorithm extracting p-SCC.

Digraph kernelization

In the previous part, we explained the general idea of our approach, finding communities cores by searching p-SCC, which means exploring a digraph in order to find circuits. Before describing the whole method based on core detection, we show how we can optimize the search for p-SCCs by "cleaning" the digraph, meaning excluding nodes and arcs which should never belong to a circuit: this brings us to the notion of graph kernelization. We are interested in the kernelization which is used in the FVS problem (Feedback Vertex Set) [START_REF] Thomassé | A quadratic kernel for feedback vertex[END_REF]. Its purpose is to find a set of nodes which gives a graph without cycles if this set of nodes (and their adjacent edges) is removed from the graph. The general interest of the kernelization process is the reduction of the size of the input given to an algorithm which is not polynomial in time (in most cases): the graph is being "compressed". As we work in the directed case, we used the kernelization technique applied to the directed FVS [START_REF] Fleischer | Experimental study of fpt algorithms for the directed feedback vertex set problem[END_REF], following the four first rules of the method. The two first rules imply a simple digraph (no self-loop, no multiple arcs), the third one removes isolated nodes (degree equal to zero), and the last one removes the chained nodes, by removing nodes from the digraph while it contains nodes with only one outgoing or incoming arc. The experiment part shows that the kernelization operation can be very effective on real-world networks.

Digraph cores decomposition method

This section describes our method for core detection in directed networks. Let use the following notations: G is the input digraph (network), and K is the set of output cores. The method follows these steps, considering a given p:

1. Kernelize G. 2. For each node of G as the starting node, process p-SCC. 3. Sort p-SCCs by size. Starting from the biggest one, put them one by one in K if it doesn't intersect existing cores already inserted into K. In case of cores having the same size, take the most connected one (biggest amount of arcs). 4. (optionnal) Remove p-SCCs with size inferior to a given threshold K min .

Illustration:

The figure 3 gives an illustration of our method, step by step, with p = 4 (meaning we search for circuits of length at most 3). Let take a digraph (a), and apply the first step which is kernelization (b). Some nodes are ignored, and won't be considered. The second step processes 4-SCCs, node by node: in the example (c), only five iterations are represented (nodes with labels 4 ,5 ,8 ,10 ,13), and for each node, a 4-SCC is computed, which can be the same for several nodes (nodes 5 and 8 produce the same 4-SCC, same thing for nodes 4 and 13). The last step (d) extracts the biggest and non-intersecting 4-SCCs giving the final result with 3 cores. This method returns a set of cores which can be used to cluster the rest of the network to have a complete clustering. As some clustering methods like k-means algorithm, the number of communities is set by the number of cores. In this article, as we don't focus on clustering methods, our experiments use a simple aggregative method like center-based clustering methods, and assign each node to the community having the nearest core. The distance corresponds to the chain length from node to core (at this point, link orientation is ignored). At each step, the cores absorb their nearest neighbors, and the process continues until all the nodes belong to a group. If a node can be absorbed by more than one community, meaning the distance between several cores is equal, it will be absorbed by the community sharing the biggest amount of edges with it. The K min value can be useful to avoid too small cores that shouldn't be considered.

Experiments

Validating communities structures corresponds to validate a clustering method. The difficulty is to find an objective measure of quality of clusters. For our experiments, we use V-Measure which is an alternative to F-Measure, and Normalized Mutual Information from the information theory field to compare the clustering obtained by our method to the reference classes. We based our experiments on real data networks, as some experiments have already been done on generated networks (LFR Benchmarks [START_REF] Lancichinetti | Benchmarks for testing community detection algorithms on directed and weighted graphs with overlapping communities[END_REF]) in our previous work with good results [START_REF] Levorato | Detection of communities in directed networks based on strongly p-connected components[END_REF]. Moreover, there is an issue with the LFR Benchmark, as it produces graphs already kernelized, which puts a strong constraint on generated graphs. On the contrary, in the experiment part, we observe that the real-world networks we used are strongly kernelizable.

Clustering Evaluation

The entropy notion is used to express the used measures, and is noted as follow, with X and Y two discrete random variables: H(X) and H(Y ) for the marginal entropies, H(X|Y ) and H(Y |X) for the conditional entropies, H(X, Y ) for the joint entropy. Measures give results between 0 (worst matching) and 1 (best matching). Here are the two evaluation measures definitions:

-V-Measure [START_REF] Rosenberg | V-measure: A conditional entropy-based external cluster evaluation measure[END_REF] is an entropy based-evaluation measure, composed of two concepts: completeness and homogeneity. The V-Measure corresponds to the harmonic mean of these two concepts, just like F-Measure with precision and recall. With C a set of classes (reference), K a set of clusters (unsupervised method), the homogeneity is defined as:

h = 1 if H(C, K) = 0 1 -H(C|K) H(C) else
and the completeness is defined as:

c = 1 if H(K, C) = 0 1 -H(K|C) H(K) else
A clustering result satisfies homogeneity if all of its clusters contain only elements which are members of a single class, and a clustering result satisfies completeness if all the elements that are members of a given class are elements of the same cluster. The V-Measure is based on homogeneity and completeness scores such as:

V β = (1 + β) • h • c (β • h) + c)
with a β parameter which can be used to weight homogeneity and completeness scores.

This measure solves the problem of other measures like F-Measure which only consider the contributions from the clusters which match a target class, and has sub-measures useful to evaluate precisely the quality of a clustering solution. In our experiments, we set β = 1 (balanced weights).

-Normalized Mutual Information (NMI) [START_REF] Danon | Comparing community structure identification[END_REF] Mutual information is a measure used in information theory domain, giving the amount of information that one random variable contains about another. The entropy reflects the self-information of a random variable. Mutual information can be can be summarized in a measure of the distance between two probability distributions. This measure is defined between the cluster assignments K and a pre-existing labeling set of classes C normalized by:

N M I(K, C) = I(K, C) H(K)H(C)
with I(K, C) the mutual information of K and C such that:

I(K, C) = H(K) -H(K|C)

Results

In order to test our approach, we used directed network datasets already known in the literature [START_REF] Lada | The political blogosphere and the 2004 u.s. election: divided they blog[END_REF][START_REF] Sen | Collective classification in network data[END_REF]. Three networks have been used: a political blog network about US politics, and two citation networks (Cora and Citeseer) (see tab. 1). Our method have been applied on the biggest connected component (WCC) to avoid noisy results. Several characteristics of the network are given like density, degree information, communities maximum and minimum sizes, but also the mixing parameter µ [START_REF] Lancichinetti | Community detection algorithms: A comparative analysis[END_REF] and the directed modularity Q d [START_REF] Leicht | Community structure in directed networks[END_REF]. The mixing parameter is used by Lancichinetti and al. to generate datasets of networks. It represents the fraction of all links in a particular community that end outside this community. More its value increases, less the communities are well defined and less easily detectable. Network modularity is defined as the fraction of all links that lie within communities minus the expected value of the same quantity in a network in which nodes have the same degree but with links placed randomly. It gives an idea on the "good" separation of communities. For each network, we seek for cores, and study the results using only V-measure. Then we focus on the communities result with both NMI and V-Measure.

Core decomposition

The quality of the core depends on the kernelization process, as graph compression is highly different from a network to another (tab. 2). In tab. 3, we compare the obtained cores with the reference classes, using the p parameter which corresponds to the path length of a p-SCC, and the K min parameter which is the minimum core size (only relevant results are shown). In our experiments, the use of the minimum core size has been relevant. In most cases, cores have a good completeness score, meaning that we succeed in having nodes which belong to a single class in only one core. On the other hand, the homogeneity score tends to be better when the threshold of the minimum core size is increased (Political Blog and Cora networks), having nodes of a same core belonging to a single class. The interpretation that can be made from these results is that the more the graph is compressed, the less the K min gets an high value. When too many nodes are available to build cores, the K min threshold has to be high to remove some eventual noise, giving less nodes usable in the core creation (illustration fig. 4). In the results of tab. 3, we first consider completeness value, and then the homogeneity value. We give more importance to the completeness score, as it gives better results in the final process communities detection. Communities detection Using an aggregative method, the cores first absorb nodes which are in the kernel but not in the cores, giving pre-built communities. Then, the nodes outside the kernel are absorbed by the pre-built communities to give a final clustering of communities. Clusters are not strongly connected, but unilaterally connected. The results in tab. [START_REF] Cheng | Predicting reciprocity in social networks[END_REF] show that even with a naive method of clustering, the communities structures remain "acceptable". The cores used in the final clustering process are the cores having the best completeness scores in the core decomposition operation. Compared to InfoMap algorithm (tab. 5), which stays the method with best results obtained in directed networks [START_REF] Lancichinetti | Community detection algorithms: A comparative analysis[END_REF], our results look pretty similar, with best performance on the Political Blogs network. Infomap also tends to subdivide too much the networks, finding 30 to 60 times more communities than expected. Observing the results, we can make the assumption that the compression of graphs impacts the quality of cores, and therefore the detection of communities. With a small amount of nodes in the kernel, the choice to make between the nodes to build the cores is important, as it determines the final process of communities detection. Also, having only big cores means setting a too high K min threshold value, which can have a negative impact on the cores detection, and some communities cannot be found in the process. For instance, in fig. 5, the central community has been found by our method using the parameters p = 4, K min = 2. If we set K min = 3, cores of size 2 are excluded, and this central community is not detected.

Conclusion

In this article, we focused on an approach dedicated to directed networks, and we gave a method allowing the decomposition of these networks into cores. These cores can be used by any clustering method based on centers to detect communities. Our various contributions can be presented as follows:

-Using the strongly p-connected components is the key of our approach, which has not been used in any scientific work in order to detect cores or communities in directed networks. Moreover, we provide a simple and efficient algorithm to generate these p-SCCs in a digraph. This approach can be classified in the pattern identification category that we can find in some method classification, while being flexible enough. -The interest of using kernelization process has been highlighted : it reduces the core detection process, and can give some information on the network structure. The hypothesis of having the compression rate of a digraph correlated to the mixing parameter has to be taken into account, and more experiments on real-world networks have to be done to verify this intuition. -During our experiments, we obtained some encouraging results on the network datasets we used. An important thing about these results is that we didn't take into account the modularity concept in our approach. As a large part of the communities detection algorithms are dedicated to modularity optimization [START_REF] Labatut | Detection and interpretation of communities in complex networks: Practical methods and application[END_REF], we want to stress the point that we can have interesting results in communities detection without this concept, which some limits are already known even in the undirected case [START_REF] Lancichinetti | Limits of modularity maximization in community detection[END_REF].

Several options can be considered for the continuation of this work. As we said, we have to apply our method to others real-world datasets. We should also study how to increase the quality of the core detection, and it could be interesting to have the possibility to automatically fix the K min threshold value. Testing other based-centers clustering methods should be done too. Also, the case of overlapping communities should be considered, as our approach could be quickly adaptable with p-SCCs which naturally overlap each other. In our opinion, our work points out that no clear or unanimous consensus about the definition of communities exists, and provides a new point of view on the detection of communities into directed networks, being omnipresent in the Web nowadays.

Fig. 1 .

 1 Fig. 1. Examples of p-SCCs: (a): nodes are connected by paths of length at most 2 (2-SCC) (b): nodes are connected by paths of length at most 3 (3-SCC).

Fig. 3 .

 3 Fig. 3. Illustration of the decomposition method of a digraph into cores.

Fig. 4 .

 4 Fig. 4. Illustration of core decomposition in the Political Blog network.

  (a) Clustering using an aggregative method based on core decomposition. (b) Reference classes.

Fig. 5 .

 5 Fig. 5. Communities detection comparison on the Cora citation network (p = 4, Kmin = 2).

  Input: G: digraph, s:starting node, p:path length (even integer) Data: astack: stack of arcs, vpath: stack of nodes, c: integer (circuit size) Remark: Aout(k) represents set of outing arcs of the node k. source(a) represents the source node of the arc a. dest(a) represents the destination node of the arc a.

	end
	vpath.push(w);
	end
	z ← dest(a);
	if z = s then
	C ← C ∪ vpath;
	end
	else
	if |vpath| < c then
	foreach b in Aout(z) do
	astack.push(b);
	end
	vpath.push(z);
	end
	end
	end
	return C;

Result: C:set of nodes (p-SCC) C ← ∅; C ← C ∪ {s}; vpath.push(s); c ← (p+2) 2 ; foreach a in Aout(s) do astack.push(a); end while astack = ∅ do a ← astack.pop(); w ← vpath.peek(); if source(a) = w then while source(a) = w do w ← vpath.pop();

Table 1 .

 1 Network datasets.

	Directed Network Political Blog	Cora	Citeseer
	|V |	1,222	2,485	2,120
	|A|	19,024	5,209	3,768
	Classes	2	7	6
	Density	1.27%	0.08%	0.08%
	Degrees	kmean = 31 kmean = 4	kmean = 4
		kmin = 1	kmin = 1	kmin = 1
		kmax = 467 kmax = 169 kmax = 100
	Communities size |C|min = 588 |C|min = 131 |C|min = 115
		|C|max = 636 |C|max = 726 |C|max = 532
	µ	0.09	0.18	0.28
	Q d	0.41	0.63	0.51

Table 2 .

 2 Network kernel sizes.

	Directed Network	Political Blog Cora Citeseer
	|V |K	811	399	69
	|A|K	15,833	786	97
	Compression rate (nodes)	33%	84% 97%

Table 3 .

 3 Cores detection on real-world networks.

	Network			Measures	Amount of
					Communities
		h	c	V	NMI
	Political Blog 0.70385 0.69929 0.70156 0.70116	2
	Cora	0.35335 0.46349 0.40099 0.40469	28
	Citeseer	0.28162 0.38734 0.32613 0.32742	26

Table 4 .

 4 Real-world networks communities detection results based on core decomposition.

	Network		Measures		Amount of
						Communities
		h	c	V	NMI
	Political Blog 0.281	0.759	0.410	0.462	56
	Cora	0.259	0.833	0.395	0.465	438
	Citeseer	0.228	0.721	0.347	0.404	311

Table 5 .

 5 Real-world networks communities detection results obtained by InfoMap.

Acknowledgments

We would like to thank Pr. Ioan Todinca, Mathieu Liedloff, Anthony Perez and the GAMoC team of the LIFO at University of Orléans (FR) for their listening and advices. Special thanks to Guillaume Cleuziou from the Machine Learning team.