B Galerne
email: galerne@cmla.ens-cachan.fr

Y Gousseau
email: gousseau@telecom-paristech.fr

The Transparent Dead Leaves Model

Keywords: Germ-grain model, dead leaves model, transparency, occlusion, image modeling, texture modeling

This paper introduces the transparent dead leaves (TDL) random field, a new germ-grain model in which the grains are combined according to a transparency principle. Informally, this model may be seen as the superposition of infinitely many semi-transparent objects. It is therefore of interest in view of the modeling of natural images. Properties of this new model are established and a simulation algorithm is proposed. The main contribution of the paper is to establish a central limit theorem, showing that when varying the transparency of the grain from opacity to total transparency, the TDL model ranges from the dead leaves model to a Gaussian random field.

Introduction

This paper deals with the stochastic modeling of physical transparency. The main contribution is the introduction and study of a new germ-grain model in which the grains are combined according to a transparency principle. To the best of our knowledge, this type of interaction between grains has not been studied before. Classical interactions between grains include addition for shot-noise processes [START_REF] Rice | On generalized shot noise[END_REF][START_REF] Heinrich | Normal convergence of multidimensional shot noise and rates of this convergence[END_REF], union for Boolean models [START_REF] Stoyan | Stochastic geometry and its applications second[END_REF][START_REF] Schmitt | Estimation of the density in a stationary Boolean model[END_REF], occlusion for dead leaves models [START_REF] Matheron | Schéma booléen séquentiel de partition aléatoire[END_REF][START_REF] Jeulin | Dead leaves models: from space tesselation to random functions[END_REF][START_REF] Bordenave | The dead leaves model: a general tessellation modeling occlusion[END_REF] or multiplication for compound Poisson cascades [START_REF] Barral | Multifractal products of cylindrical pulses[END_REF][START_REF] Chainais | Infinitely divisible cascades to model the statistics of natural images[END_REF].

The proposed model, that we call transparent dead leaves (TDL), is obtained from a collection of grains (random closed sets) indexed by time, as for the dead leaves model of G. Matheron. We assume that each grain is given a random gray level (intensity). Informally, the TDL model may be seen as the superposition of transparent objects associated with the grains. When adding a new grain, new values are obtained as a linear combination of former values and the intensity of the added grain, as illustrated in Figure 1. More precisely, the superposition of a transparent grain X with gray level a on an image (a function f : R d → R) results in a new image f defined for each y ∈ R d by: where α ∈ (0, 1] is a transparency coefficient. The TDL model is defined as the sequential superposition of grains of a marked Poisson point process i δ (ti,xi,Xi,ai) , with i δ (ti,xi) a homogeneous Poisson point process in (-∞, 0) × R d and X i , a i , i.i.d. random sets and random variables respectively. In particular, the value of the TDL at each point results from the superposition of infinitely many semi-transparent objects.

f (y) = αa + (1 -α)f (y) if y ∈ X, f (y) otherwise, (1)
The main motivation to define such a model originates from the modeling of image formation. Indeed, natural images are obtained from the light emitted by physical objects interacting in various ways. In the case of opaque objects, the main interaction is occlusion. That is, objects hide themselves depending on their respective positions with respect to the eye or the camera. A simple stochastic model for occlusion is given by the dead leaves model, which is therefore useful for the modeling of natural images [START_REF] Gousseau | Modeling occlusion and scaling in natural images[END_REF][START_REF] Cao | Dead leaves model for measuring texture quality on a digital camera[END_REF]. When objects are transparent, their interaction may be modeled by Formula (1). This is well known in the field of computer graphics, see [START_REF] Enderton | Stochastic transparency[END_REF] where the same principle is used for the creation of synthetic scenes. In this case, transparency is a source of heavy computations, especially in cases where objects are numerous (typically of the order of several thousands), e.g. in the case of grass, fur, smoke, fabrics, etc. The transparency phenomenon may also be encountered in other imaging modality where images are obtained through successive reflexion-transmission steps, as in microscopy or ultrasonic imaging. A related non-linear image formation principle is at work in the field of radiography. In such cases, it is useful to rely on accurate stochastic texture models in order to be able to detect abnormal images. The TDL may be an interesting alternative to Gaussian fields that are traditionally used, see e.g. [START_REF] Grosjean | A-contrario detectability of spots in textured backgrounds[END_REF][START_REF] Richard | Statistical tests of anisotropy for fractional Brownian textures. Application to full-field digital mammography[END_REF]. A last motivation for the TDL is that, as explained in the next paragraph, it is intermediate between the dead leaves model and Gaussian fields, two models that have proven useful for the modeling of natural textures [START_REF] Cao | Dead leaves model for measuring texture quality on a digital camera[END_REF][START_REF] Galerne | Random phase textures: Theory and synthesis[END_REF].

In this paper, we first define the TDL in Section 2 and give some elementary properties in Section 3, where we also address the problem of simulating the model and show some realizations. The TDL covariance is then computed in Section 4. Eventually, the main result of the paper is stated and proved in Section 5, namely that the normalized TDLs converge, as the transparency coefficient α tends to zero, to a Gaussian random field having the same covariance function as the shot noise associated with the grain X and with intensity one. Thus the TDLs with varying transparency coefficient α provide us with a family of models ranging from the dead leaves model to Gaussian fields.

Definition of the TDL model

As explained in the introduction, the TDL model is obtained as the superposition of transparent shapes. Formally it is defined from a marked Poisson point process, in a way similar to the dead leaves model [START_REF] Bordenave | The dead leaves model: a general tessellation modeling occlusion[END_REF]. Let F denote the set of closed subsets of R d . On the state space

S = (-∞, 0) × R d × F × R,
equipped with its natural product σ-algebra, we define the point process

Φ = i δ (ti,xi,Xi,ai) , (2)
where

• {(t i , x i )} is a stationary Poisson point process of intensity 1 in the half space (-∞, 0) × R d ,
• (X i ) i is a sequence of i.i.d. random closed sets (RACS) with distribution P X which is independent of the other random objects,

• (a i ) i is a sequence of i.i.d. real random variables with distribution P a which is also independent of the other random objects.

Equivalently, Φ is a Poisson point process with intensity measure µ = λ⊗ν d ⊗P X ⊗P a , where λ denotes the restriction of the one-dimensional Lebesgue measure to (-∞, 0) and ν d denotes the d-dimensional Lebesgue measure on R d . Each point (t i , x i , X i , a i ) of the Poisson process Φ is called a leaf. Having fixed a transparency coefficient α ∈ (0, 1], the TDL process f is obtained by sequentially combining the elements of Φ according to Formula (1), which results in the following definition.

Definition 1 (Transparent Dead Leaves model). The Transparent Dead Leaves model with transparency coefficient α associated with the Poisson process Φ defined by Equation (2) is the random field f : R d → R defined by

f (y) = i∈N 1 (y ∈ x i + X i ) αa i (1 -α) ( j∈N 1(tj∈(ti,0) and y∈xj +Xj )) . (3)
Let us justify that Formula (3) agrees with the informal description of the TDL model. According to Equation (1), the impact of the leaf (t i , x i , X i , a i ) is to add αa i and to attenuate the previous contributions by a factor (1-α). Hence, the contribution of the leaf (t i , x i , X i , a i ) at a point point y ∈ x i + X i is αa i multiplied by (1 -α) to the number of leaves fallen on the point y after the leaf (t i , x i , X i , a i ), that is after time t = t i . This number is precisely the exponent of (1 -α) in Formula (3): j∈N 1 (t j ∈ (t i , 0) and y ∈ x j + X j ) .

Remark (Random functional). Denoting by N(S) the set of point processes taking value in the state space S, one remarks that the TDL random field f (y) has the form

f (y) = (ti,xi,Xi,ai)∈Φ g(y, (t i , x i , X i , a i ), Φ),
where g : R d × S × N(S) → R is a measurable function given by Formula (3). Similar random functionals interpreted as the sum of contributions from each point of a (possibly marked) point process {x i } ⊂ R d appear in several contexts in stochastic geometry. In particular, general central limit theorems hold when the intensity of the point process {x i } tends to infinity, see e.g [START_REF] Baryshnikov | Gaussian limits for random measures in geometric probability[END_REF][START_REF] Penrose | Gaussian limits for random geometric measures[END_REF]. Note however that our framework is different since the Poisson process has an additional time component t i , and consequently, there is always an infinite number of leaves (t i , x i , X i , a i ) influencing the value f (y) (as will be clarified by Proposition 3).

Remark (Variable transparency). For the sake of simplicity, the transparency parameter α is assumed to be the same for all objects. However, one may attach a random transparency α i to every object in Definition 1 and generalize the results of Sections 3 and 4, as will be briefly commented thereafter.

Since the distribution of the Poisson process Φ is invariant under shifts of the form (t, x, X, a) → (t, x + y, X, a), the TDL f is a strictly stationary random field.

Before establishing further properties of the TDL random field f , let us introduce some notation and specify several assumptions.

Notation: One defines β = 1 -α and one denotes by X and a, respectively, a RACS with distribution P X and a r.v. with distribution P a which are both independent of all the other random objects. In addition, γ X denotes the mean geometric covariogram of the RACS X, that is the function defined by γ X (τ ) = E (ν d (X ∩ (τ + X))), τ ∈ R d (we refer to [START_REF] Matheron | Random sets and integral geometry[END_REF][START_REF] Lantuéjoul | Geostatistical simulation: models and algorithms[END_REF] for properties of the mean geometric covariogram).

Assumptions: Throughout the paper, it is assumed that

0 < E (ν d (X)) < +∞.
This hypothesis ensures that each point y ∈ R d is covered by a countable infinite number of leaves of Φ, whereas the number of leaves falling on y during a finite time interval [s 1 , s 2 ] is a.s. finite. We also assume that E(a 2 ) < +∞.

3 One-dimensional marginal distribution and simulation of the TDL model

The Poisson process of the leaves intersecting a set

As one can observe from Equation (3), the only leaves which have a contribution to the sum defining f (y) are the leaves (t i , x i , X i , a i ) such that y ∈ x i +X i . When considering the restriction of f to a Borel set G the only leaves of interest are the ones intersecting G, i.e. the leaves (t i , x i , X i , a i ) such that x i + X i ∩ G = ∅. The next proposition characterizes the distribution of such leaves, a result to be used further in the paper. We first introduce two notation: if A and B are two Borel sets then Ǎ = {-x : x ∈ A} and A ⊕ B = {x + y : x ∈ A and y ∈ B}.

Remark that x + X ∩ G = ∅ ⇐⇒ x ∈ G ⊕ X.
Proposition 2 (The Poisson process of the leaves intersecting a Borel set). Let G ⊂ R d be a Borel set such that 0 < E ν d X ⊕ Ǧ < +∞ and let Φ be the Poisson process on S = (-∞, 0) × R d × F × R with intensity measure µ = λ ⊗ ν d ⊗ P X ⊗ P a . Denote by Φ G the point process of the leaves of Φ which intersect G, that is

Φ G = {(t, x, X, a) ∈ Φ : x + X ∩ G = ∅} ,
and denote by

A G ⊂ R d × F the set A G = {(x, X) : x + X ∩ G = ∅}. Then Φ G is a
Poisson process on S with intensity measure

µ G = λ ⊗ (ν d ⊗ P X ) A G ⊗ P a .
It is an independently marked Poisson process with ground process

Π G = t : (t, x, X, a) ∈ Φ G , a homogeneous Poisson process on (-∞, 0) of intensity E ν d X ⊕ Ǧ , and with mark distribution 1 E ν d X ⊕ Ǧ (ν d ⊗ P X ) A G ⊗ P a .
Proof. Φ G is the restriction of the Poisson process Φ to the measurable set

(t, x, X, a) ∈ (-∞, 0) × R d × F × R : (x, X) ∈ A G ,
thus Φ G is a Poisson process and its intensity measure µ G is the restriction of µ to the above set. As for the interpretation of Φ G as an independently marked one-dimensional Poisson process, it is based on the factorization of the intensity measure µ G (see [START_REF] Baddeley | Spatial point processes and their applications[END_REF]Section 1.8] or [START_REF] Schneider | Stochastic and Integral Geometry[END_REF]Section 3.5]). Indeed we have

0 < ν d ⊗ P X A G = F R d 1 y ∈ G ⊕ Y ν d (dy)P X (dY ) = E ν d X ⊕ Ǧ < +∞,
and thus we can write

µ G = E ν d X ⊕ Ǧ λ ⊗ 1 E ν d X ⊕ Ǧ (ν d ⊗ P X ) A G ⊗ P a ,
where the measure between square brackets is a probability distribution.

One-dimensional marginal distribution

Proposition 3 (One-dimensional marginal distribution). Let y be a point in R d . Then there exists a sequence (a(y, k)) k∈N of i.i.d. r.v. with distribution P a such that

f (y) = α +∞ k=0 a(y, k)β k .
In particular we have

E(f (y)) = E(a) and Var (f (y)) = α 2 -α Var(a).
Informally, a(y, k) is the color of the (k + 1)-th leaf falling on y. (a(y, k)) k∈N is thus an ordered subfamily of the r.v. (a i ) i∈N which depends on y.

Proof. According to Proposition 2 the point process Φ {y} of the leaves which cover y is an independently marked Poisson process, the ground process of which is a Poisson process on (-∞, 0) with intensity 0 < E(ν d (X)) < +∞. Hence the falling times of the leaves of Φ {y} are a.s. distinct and we can number (in a measurable way [START_REF] Schneider | Stochastic and Integral Geometry[END_REF]) the leaves (t(y, k), x(y, k), X(y, k), a(y, k)), k ∈ N, according to an anti-chronological order:

0 > t(y, 0) > t(y, 1) > t(y, 2) > . . .
Proposition 2 also gives the distribution of the marks (x(y, k), X(y, k), a(y, k)), and in particular it shows that the r.v. a(y, k), k ∈ N are i.i.d. with distribution P a . As already mentioned, the only leaves involved in the sum which defines f (y) are the leaves of Φ {y} . Besides, using the above numbering we have for all k ∈ N (tj ,xj ,Xj ,aj )∈Φ

1 (t j ∈ (t(y, k), 0) and y ∈ x j + X j ) = k.
Hence Equation (3) becomes

f (y) = α +∞ k=0 a(y, k)β k ,
and the result follows.

Remark (Influence of the transparency coefficient α). Let us write f α for the TDL model with transparency coefficient α ∈ (0, 1]. Proposition 3 shows that the expectation of f α does not depend on α. In contrast, the variance Var (f α (y

)) = α 2 -α Var(a)
decreases as α decreases. Besides, Var (f α (y)) tends to 0 as α tends to 0 (recall that the model is not defined for α = 0). However, a central limit theorem for random geometric series [START_REF] Bovier | A law of the iterated logarithm for random geometric series[END_REF] shows that for all y ∈ R d the family of r.v.

f α (y) -E (f α ) Var (f α )
α converges in distribution to a standard normal distribution as α tends to 0. This pointwise convergence result will be extended in Section 5, where it will be shown that the family of normalized random fields y →

f α (y) -E (f α ) Var (f α ) α converges in the sense of finite-dimensional distributions.

Simulation of the TDL model

In this section we draw on Proposition 3 to obtain a simulation algorithm for the restriction of the TDL model f to a finite set U ⊂ R d (e.g. a finite grid of pixels). The algorithm is based on a coupling from the past procedure, as the algorithm developed by Kendall and Thönnes [START_REF] Kendall | Perfect simulation in stochastic geometry[END_REF] for simulating the dead leaves model (see also [START_REF] Jeulin | Dead leaves models: from space tesselation to random functions[END_REF][START_REF] Lantuéjoul | Geostatistical simulation: models and algorithms[END_REF]). This algorithm consists in sequentially superimposing transparent random objects but, contrary to the forward procedure described by Equation (1), each new object is placed below the former objects. In the case of the dead leaves model, this yields a perfect simulation algorithm. For the TDL model f , simulation is not perfect since the values f (y) are the limits of convergent series. Nevertheless, supposing that the intensities a i are bounded, we propose for any precision ε > 0 an algorithm which produces an approximation f of f . This approximation satisfies

P sup y∈U f (y) -f (y) ≤ ε = 1
therefore providing a kind of perfect simulation with precision ε > 0.

In the remaining of this section we suppose that the colors a i are a.s. bounded by A > 0. The control of the precision is based on the following elementary lemma. Lemma 4 shows that to approximate f (y) with a tolerance ε > 0 it is enough to cover the point y with (at least) N (ε) leaves, where

N (ε) is the smallest integer n such that Aβ n ≤ ε, that is N (ε) = log (ε/A) log(β)
. The following simulation algorithm relies on this observation.

Algorithm 1 (Simulation of the TDL model with tolerance ε > 0). Let U ⊂ R d be a finite set. Given a precision ε > 0, an approximation f of the TDL model f is computed by controlling the number of leaves L at each point:

• Initialization: For all y ∈ U , f (y) ← 0; L(y) ← 0;

• Computation of the required number of leaves: N (ε) = log (ε/A) log(β) ;

• Iteration: While inf y∈U L(y) < N (ε) add a new leaf:
1. Draw a leaf (x, X, a) hitting U :

(a) Draw X ∼ P X ; (b) Draw x uniformly in U ⊕ X; (c) Draw a ∼ P a ; 2. Add the leaf (x, X, a) to f : for all y ∈ U , f (y) ← f (y)+1 (y ∈ x + X) αaβ L(y) ; 3. Update the leaves layer L: for all y ∈ U , L(y) ← L(y) + 1 (y ∈ x + X) ;
Clearly Algorithm 1 a.s. terminates if every point of U is covered by N (ε) leaves in an a.s. finite time. This is always the case since U is a finite set and E (ν d (X)) > 0.

Several realizations of some TDL models are represented in Figure 2. Remark that as soon as α < 1, the TDL random field is not piecewise constant: any region is intersected by the boundaries of some leaves, producing discontinuities.

Covariance of the TDL model

This section is devoted to the computation of the covariance of the TDL. A classical way to achieve this would be to use Palm calculus, leading to relatively heavy computations in this case. Instead, we chose an alternative way relying on some memoryless property of the TDL, as explained below.

The following proposition is an extension of the fact that if 0 > t 0 > t 1 > t 2 > . . . is a homogeneous Poisson process on (-∞, 0) then the shifted process 0 > t 1 -t 0 > t 2 -t 0 > t 3 -t 0 > . . . is also a Poisson process with the same distribution [START_REF] Kingman | Poisson Processes. Oxford Studies in Probability[END_REF]Chapter 4].

Proposition 5 (Last hitting leaf and the Poisson process preceding the last hit). Let F be a locally compact topological space with a countable base [START_REF] Schneider | Stochastic and Integral Geometry[END_REF]. Let Ψ be a Poisson process in (-∞, 0) × F with intensity measure of the form λ ⊗ η where λ is the onedimensional Lebesgue measure on (-∞, 0) and η is a measure on F . Let A ⊂ F be an η-measurable set satisfying 0 < η(A) < +∞. Define

t 0 = sup {t i |(t i , y i ) ∈ Ψ ∩ ((-∞, 0) × A)} ,
y 0 the a.s. unique y ∈ F such that (t 0 , y) ∈ Ψ ∩ ((-∞, 0) × A), and

Ψ t0 = (ti,yi)∈Ψ 1 (t i < t 0 ) δ (ti-t0,yi) .

Then

• t 0 , y 0 , and Ψ t0 are mutually independent,

• -t 0 has an exponential distribution with parameter η(A),

• y 0 has distribution Q A defined for all B ∈ B(F ) by Q A (B) = η(B ∩ A) η(A) ,
• Ψ t0 is a Poisson process with intensity measure λ ⊗ η, i.e. Ψ t0 has the same distribution as Ψ.

Proposition 5 will be applied below to the Poisson process Φ of the colored leaves to compute some statistics of the TDL model f . As a first example, let us reobtain the expectation of f by using Proposition 5. Let y ∈ R d and let us denote by (t 0 , x 0 , X 0 , a 0 ) the leaf which hits y at the maximal time t 0 . Then one can decompose f (y) into

f (y) = αa 0 + βf t0 (y), (4)
where f t0 is the TDL model associated with the time-shifted point process Φ t0 and, as before, β = 1 -α. According to Proposition 5, a 0 has distribution P a and both point processes Φ and Φ t0 have the same distribution. Consequently, f (y) and f t0 (y) also have the same distribution, and in particular the same expectation. Hence, the above decomposition of f (y) leads to the equation

E (f (y)) = αE(a) + βE (f (y)) ,
which gives E (f (y)) = E(a), in accordance with Proposition 3.

The very same method is used below to compute the covariance of f . This method will also be applied in Section 5.3 to derive a technical result useful for the central limit theorem of Section 5.

We recall that γ X (τ ) = E (ν d (X ∩ (τ + X))) is the mean covariogram of X. In addition, one denotes the covariance of f by

Cov(f )(τ ) = Cov (f (y), f (y + τ )) = E ((f (y) -E(a)) (f (y + τ ) -E(a))) .
Proposition 6 (Covariance of the TDL model). The TDL model f is a squareintegrable stationary random field and its covariance is given by

Cov(f )(τ ) = αγ X (τ ) 2E(ν d (X)) -αγ X (τ ) Var(a), τ ∈ R d .
Proof. Let y and z be such that z -y = τ and write m = E(a) = E(f ) as a shorthand notation. We have to compute Cov(f

)(τ ) = E ((f (y) -m) (f (z) -m)) .
Denote by (t 0 , x 0 , X 0 , a 0 ) the last leaf which hits y or z at the maximal time t 0 , and let Φ t0 be the corresponding time-shifted Poisson process. According to Proposition 5, (x 0 , X 0 , a 0 ) is independent of Φ t0 . In addition Φ t0 d = Φ, and consequently, noting f t0 the TDL associated with Φ t0 , (f t0 (y), f t0 (z)) d = (f (y), f (z)). Proposition 5 also shows that a 0 has distribution P a . As for the distribution of (x 0 , X 0 ), a straightforward computation shows that

ν d ⊗ P X ({(x, X), {y, z} ∩ x + X = ∅}) = E (ν d (X ⊕ {-y, -z})) = 2γ X (0) -γ X (τ ) and ν d ⊗ P X ({(x, X), {y, z} ⊂ x + X}) = E (ν d (-y + X ∩ -z + X)) = γ X (τ ).
Hence we have

P ({y, z} ⊂ x 0 + X 0 ) = ν d ⊗ P X ({(x, X), {y, z} ⊂ x + X}) ν d ⊗ P X ({(x, X), {y, z} ∩ x + X = ∅}) = γ X (τ ) 2γ X (0) -γ X (τ ) , (5)
and by symmetry

P (y ∈ x 0 + X 0 and z / ∈ x 0 + X 0 ) = P (z ∈ x 0 + X 0 and y / ∈ x 0 + X 0 ) = γ X (0) -γ X (τ ) 2γ X (0) -γ X (τ ) .
In conditioning with respect to the coverage of the last leaf (t 0 , x 0 , X 0 , a 0 ) we have

E ((f (y) -m) (f (z) -m)) = E ((f (y) -m) (f (z) -m) |{y, z} ⊂ x 0 + X 0 ) γ X (τ ) 2γ X (0) -γ X (τ ) + E ((f (y) -m) (f (z) -m) |y ∈ x 0 + X 0 and z / ∈ x 0 + X 0 ) γ X (0) -γ X (τ ) 2γ X (0) -γ X (τ ) + E ((f (y) -m) (f (z) -m) |z ∈ x 0 + X 0 and y / ∈ x 0 + X 0 ) γ X (0) -γ X (τ ) 2γ X (0) -γ X (τ ) .
By symmetry it is clear that the two last terms of the above sum are equal. On the event {{y, z} ⊂ x 0 + X 0 } we have

f (y) -m = α(a 0 -m) + β (f t0 (y) -m) and f (z) -m = α(a 0 -m) + β (f t0 (z) -m) , so that (f (y) -m) (f (z) -m) = α 2 (a 0 -m) 2 + β 2 (f t0 (y) -m) (f t0 (z) -m) + αβ(a 0 -m) ((f t0 (y) -m) + (f t0 (z) -m)) .
By Proposition 5, a 0 , (x 0 , X 0 ), and (f t0 (y), f t0 (z)) are mutually independent, hence

E ((f (y) -m) (f (z) -m) |{y, z} ⊂ x 0 + X 0 ) = α 2 E (a 0 -m) 2 + β 2 E ((f t0 (y) -m) (f t0 (z) -m)) = α 2 Var(a) + β 2 Cov (f (y), f (z)) .
On the event {y ∈ x 0 + X 0 and z / ∈ x 0 + X 0 } we have

f (y) -m = α(a 0 -m) + β (f t0 (y) -m) and f (z) -m = f t0 (z) -m.
Hence, by the same arguments,

E ((f (y) -m) (f (z) -m) |y ∈ x 0 + X 0 and z / ∈ x 0 + X 0 ) = β Cov (f (y), f (z)) .
Replacing the terms in the decomposition of E ((f (y) -m) (f (z) -m)) leads to an equation involving the covariance Cov (f (y), f (z)), the values γ X (0) and γ X (τ ) of the mean covariogram of X, and the variance Var(a). Simplifying this equation one obtains the enunciated formula.

Remark (Variable transparency and second-order property). The technique used in this section enables us to generalize second-order formulas to the case where the transparency parameter α is assumed to be different for each object, that is, when it is assumed that each object X i is assigned a transparency α i distributed as a random variable α and independent of other objects. First, it is straightforward to show that in this case we still have E(f (y)) = E(a). Then, a simple application of Formula (4) yields Var f (y) = E(α 2 ) Var(a)(2E(α) -E(α 2 )) -1 . Observe that a direct computation starting from the definition of f would be much more arduous. Eventually, applying the same technique, one can show that the covariance of the model with variable transparency satisfies, for

τ ∈ R d , Cov(f )(τ ) = E(α 2 )γ X (τ ) 2E(α)E(ν d (X)) -E(α 2 )γ X (τ )
Var(a).

Gaussian convergence as the objects tend to be fully transparent

Recall that the TDL model with transparency coefficient α is denoted f α .

Theorem 7 (Normal convergence of the TDL model). Suppose that Var(a) > 0.

Then, as the transparency coefficient α tends to zero, the family of random fields

f α -E (f α )
Var (f α ) α converges in the sense of finite-dimensional distributions to a stationary Gaussian random field with covariance function

C(τ ) = γ X (τ ) E(ν d (X)) = γ X (τ ) γ X (0) .
Before proving Theorem 7 let us illustrate the normal convergence of the normalized family of r.v.

f α -E (f α )
Var (f α ) α with Figure 3. The five first images of Figure 3 are normalized TDL realizations obtained from the same random colored leaves but with various transparency coefficients α. The last image is a realization of the limit Gaussian random field given by Theorem 7. Observe that this Gaussian field is also the limit of the normalized shot noise associated with X when the intensity of germs tends to infinity [START_REF] Heinrich | Normal convergence of multidimensional shot noise and rates of this convergence[END_REF]. The remaining of this section is devoted to the proof of Theorem 7. The proof consists in showing that the finite moments of the normalized TDL random fields converge to the corresponding moments of the limit Gaussian r.f. As for the computation of the covariance (see Section 4), this convergence is established by conditioning with respect to the coverage of the last leaf hitting the considered set of points (see below for details).

Some classical results of probability theory

This section gathers two classical theoretical results needed to prove Theorem 7.

Moments and convergence in distribution

Proposition 8 (Moments and convergence in distribution). Let (f n ) be a sequence of r.f. having finite moments of all order and let f G be a Gaussian r.f. If for all p ∈ N, As α decreases to 0 the normalized TDL realizations look more and more similar to the Gaussian texture 3(f).

for all (not necessarily distinct) y 1 , . . . , y p ∈ R d ,

lim n→+∞ E   p j=1 f n (y j )   = E   p j=1 f G (y j )   ,
then (f n ) converges to f G in the sense of finite-dimensional distributions.

A Recurrence relation for the moments of a multivariate normal distribution

Explicit expressions for the moments of a multivariate normal distribution are given by Isserlis' theorem that is recalled below (see e.g. [START_REF] Michalowicz | An Isserlis' theorem for mixed Gaussian variables: Application to the auto-bispectral density[END_REF] and the references therein).

Theorem 9 (Isserlis' theorem). Let Y 1 , . . . , Y 2N +1 , N ≥ 1, be normalized (i.e. E(Y i ) = 0 and Var(Y i ) = E(Y 2 i ) = 1), jointly Gaussian r.v. Then E(Y 1 Y 2 . . . Y 2N ) = E(Y i Y j ) = Cov(Y i , Y j ), and
E(Y 1 Y 2 . . . Y 2N +1 ) = 0
, where the notation means summation over all distinct ways of partitioning the set

{Y 1 , . . . , Y 2N } into N pairs {Y i , Y j } and taking the product of the N terms E(Y i Y j ) = Cov(Y i , Y j ).
From Isserlis' theorem one deduces a recurrence relation for the moments of a multivariate normal distribution. Proof. If p ≥ 2 is odd, then by Isserlis' theorem the above formula is trivial. Hence, in the following we suppose that p is even. First let j ∈ {1, . . . , p}. Factorizing with all the pairs containing j in Isserlis' identity, one obtains

E   p j=1 Y j   = p k=1 k =j Cov(Y j , Y k )E   l∈{1,...,p}\{j,k} Y l   .
The above identity is valid for all j ∈ {1, . . . , p}. Summing these p identities gives

E   p j=1 Y j   = 1 p p j=1 p k=1 k =j Cov(Y j , Y k )E   l∈{1,...,p}\{j,k} Y l   .
Now remark that in this double sum over (j, k), the terms Cov(Y j , Y k )E ( l Y l ) only depend on the pair {j, k} but not on the order. Hence, the above expression simplifies to

E   p j=1 Y j   = 2 p {j,k}⊂{1,...,p} Cov(Y j , Y k )E   l∈{1,...,p}\{j,k} Y l   .

Notation and plan of the proof of Theorem 7

Let s be a real number such that 0 < s < 1 6 (this choice for s will become clear later). For all α ∈ (0, 1], one defines the truncation operator

T α (b) =      b if b ∈ [-α -s , α -s ], α -s if b > α -s , -α -s if b < -α -s .
For all α ∈ (0, 1], f α denotes the TDL model with transparency coefficient α and

g α (y) = f α (y) -E(a)
Var(f α ) denotes its normalization. For all α ∈ (0, 1], f T α denotes the TDL model with transparency coefficient α associated with the Poisson process

Φ T = {(t i , x i , X i , T α (a i )), (t i , x i , X i , a i ) ∈ Φ},
that is the TDL model obtained by truncating the colors a i of the leaves of Φ. We have

E(f T α ) = E(T α (a)) and Var(f T α ) = α 2 -α Var(T α (a)).
As for the TDL f α , one defines

g T α (y) = f T α (y) -E(T α (a)) Var(f T α )
.

Thanks to the truncation, f T α is bounded by α -s . In particular, for all α ∈ (0, 1], f T α and g T α have finite moments of all order. We will denote by f G a centered stationary Gaussian random field with covariance function C : τ → γ X (τ ) γ X (0) . The proof of Theorem 7 is divided into two parts:

1. One shows that the normalized TDL with truncated colors g T α converges in distribution to f G by the method of moments. More precisely the sufficient condition of Proposition 8 will be shown to be true by induction on the number of points p.

2. One shows that the family g α -g T α converges to 0 in L 2 . By Slutsky's theorem (see e.g. [START_REF] Bickel | Mathematical Statistics: Basic Ideas and Selected Topics second[END_REF]), these two properties ensure that g α converges in distribution to f G .

Normal convergence of the normalized TDL having truncated colors

With the above notation, by Proposition 8, it is enough to show the following lemma.

Lemma 11 (Convergence of Moments). For all p ∈ N, for all (not necessarily distinct)

y 1 , . . . , y p ∈ R d , lim α→0 E   p j=1 g T α (y j )   = E   p j=1 f G (y j )   .
We will show this lemma by induction on p. First note that, by definition of g T α (y j ), the statement is true for p = 0 and p = 1.

For the proof by induction we now consider an integer p ≥ 2 and p points y 1 , . . . , y p of R d , and we suppose that the convergence of moments holds for all moments of order k < p.

Decomposition of the multivariate characteristic function by condi-

tioning with respect to the coverage of the last hitting leaf

We consider the random vector

g T α (y 1 ), . . . , g T α (y p ) = f T α (y 1 ) -E(T α (a)) σ T α , . . . , f T α (y p ) -E(T α (a)) σ T α ,
where σ T α = Var(f T α ). We denote by φ α (t 1 , . . . , t p ) the multivariate characteristic function of this random vector, that is

φ α (t 1 , . . . , t p ) = E e i(t1g T α (y1)+•••+tpg T α (yp)) .
One denotes by ψ α the characteristic function of the random variable T α (a)-E(T α (a)), where a follows the color distribution P a , that is ψ α (t) = E e it(Tα(a)-E(Tα(a))) .

In addition, we introduce the shorthand notation Y for the set Y = {y 1 , . . . , y p }.

In what follows we apply Proposition 5 in considering the leaves of Φ T which hit the set Y. Hence let (t 0 , x 0 , X 0 , T α (a 0 )) denote the last leaf covering at least one point of Y, and denote by g T α,t0 the corresponding time-shifted random field. Then for all y j ∈ Y, one has the decomposition g T α (y j ) = α Tα(a0)-E(Tα(a)) σ T α + βg T α,t0 (y j ) if y j ∈ x 0 + X 0 , g α,t0 (y j ) otherwise, which can also be written as follows

g T α (y j ) = α1(y j ∈ x 0 + X 0 ) T α (a 0 ) -E(T α (a)) σ T α + β 1(yj∈x0+X0) g T α,t0 (y j ).
Besides, by Proposition 5, g T α,t0 , (x 0 , X 0 ) and a 0 are mutually independent. To obtain a decomposition of the characteristic function φ α we will condition with respect to the coverage of the last leaf x 0 + X 0 . Hence, for all subsets X ⊂ Y, X = ∅, let us denote by A X ⊂ Ω the event

A X = {(x 0 + X 0 ) ∩ Y = X } and p X = P (A X ) .
The events A X , X = ∅, form a partition of the probability space Ω, and in particular

X ⊂Y,X =∅ p X = 1.
Remark that on the event A X , the above decomposition of g T α (y j ) becomes

g T α (y j ) = α1(y j ∈ X ) T α (a 0 ) -E(T α (a)) σ T α + β 1(yj∈X ) g T α,t0 (y j ).
Hence, using the mutual independence of the different random variables,

φ α (t 1 , . . . , t p ) = E e i(t1g T α (y1)+•••+tpg T α (yp)) = X ⊂Y X =∅ E e i(t1g T α (y1)+•••+tpg T α (yp)) A X p X = X ⊂Y X =∅ ψ α   α σ T α p j=1 1(y j ∈ X )t j   φ α β 1(y1∈X ) t 1 , . . . , β 1(yp∈X ) t p p X .
(6) The next step of the proof consists in differentiating the above decomposition of the multivariate characteristic function in order to obtain a recurrence relation for the moments of g T α (y 1 ), . . . , g T α (y p ) .

A Recurrence relation for the moments of g T α

We compute below the partial derivative ∂ p φα ∂t1...∂tp (t 1 , . . . , t p ) of the characteristic function φ α to obtain an expression for the moment E p j=1 g T α (y j ) . Starting from Equation ( 6), to compute ∂ p φα ∂t1...∂tp (t 1 , . . . , t p ) one needs to differentiate with respect to each variable t j the functions of the form

F X (t 1 , . . . , t p ) = ψ α   α σ T α p j=1 1(y j ∈ X )t j   φ α β 1(y1∈X ) t 1 , . . . , β 1(yp∈X ) t p .
First let us introduce some notation. In what follows, for every subset I = {i 1 , . . . , i k } ⊂ {1, . . . , p}, we write #I = k for the cardinality of I, and

∂ k f ∂t I (t 1 , . . . , t p ) = ∂ k f ∂t i1 ∂t i1 . . . ∂t i k (t 1 , . . . , t p ).
Besides, I c denotes the complementary set of indices I c = {1, . . . , p} \ I. With these notation,

∂ p F X ∂t 1 . . . ∂t p (t 1 , . . . , t p ) = p k=0 I⊂{1,...,p} #I=k ∂ k ∂t I   ψ α   α σ T α p j=1 1(y j ∈ X )t j     ∂ p-k ∂t I c φ α β 1(y1∈X ) t 1 , . . . , β 1(yp∈X ) t p = p k=0 α σ T α k I⊂{1,...,p} #I=k i∈I 1(y i ∈ X ) ψ (k) α   α σ T α p j=1 1(y j ∈ X )t j   i∈I c β 1(yi∈X ) ∂ p-k ∂t I c φ α β 1(y1∈X ) t 1 , . . . , β 1(yp∈X ) t p .
Summing over all subsets X , one has the identity

∂ p φ α ∂t 1 . . . ∂t p (t 1 , . . . , t p ) = p k=0 α σ T α k I⊂{1,...,p} #I=k X ⊂Y X =∅ i∈I 1(y i ∈ X ) ψ (k) α   α σ T α p j=1 1(y j ∈ X )t j   i∈I c β 1(yi∈X ) ∂ p-k ∂t I c φ α β 1(y1∈X ) t 1 , . . . , β 1(yp∈X ) t p p X .
Evaluating at (t 1 , . . . , t p ) = (0, . . . , 0), this gives

∂ p φ α ∂t 1 . . . ∂t p (0, . . . , 0) = p k=0 α σ T α k I⊂{1,...,p} #I=k X ⊂Y X =∅ i∈I 1(y i ∈ X ) ψ (k) α (0) i∈I c β 1(yi∈X ) ∂ p-k ∂t I c φ α (0, . . . , 0) p X = p k=0 α σ T α k ψ (k) α (0) I⊂{1,...,p} #I=k ∂ p-k ∂t I c φ α (0, . . . , 0)     X ⊂Y X =∅ i∈I 1(y i ∈ X ) i∈I c β 1(yi∈X ) p X     .
In the above sum, remark that for k = 0, I = ∅ and thus all the terms are proportional to ∂ p φα ∂t1...∂tp (0, . . . , 0). Besides, since T α (a)-E(T α (a)) is centered, ψ

α (0) = 0, and thus for k = 1 all the terms are zero. Hence we have the following equation:

∂ p φ α ∂t 1 . . . ∂t p (0, . . . , 0)     1 - X ⊂Y X =∅   p j=1 β 1(yj∈X )   p X     (7) = p k=2 α σ T α k ψ (k) α (0) I⊂{1,...,p} #I=k ∂ p-k ∂t I c φ α (0, . . . , 0)     X ⊂Y X =∅ i∈I 1(y i ∈ X ) i∈I c β 1(yi∈X ) p X     .

Recurrence relation for the limit of the moments

The next step of the proof consists in dividing by α and letting α tend to 0 in Equation [START_REF] Cao | Dead leaves model for measuring texture quality on a digital camera[END_REF] above. First, recalling that β = 1 -α, and using that

p X = 1, one has 1 - X ⊂Y X =∅   p j=1 β 1(yj∈X )   p X = X ⊂Y X =∅ p X - X ⊂Y X =∅ β #X p X = X ⊂Y X =∅ 1 -(1 -α) #X p X . Hence lim α→0 1 α X ⊂Y X =∅ 1 -(1 -α) #X p X = X ⊂Y X =∅ (#X )p X ,
and by definition of p X ,

X ⊂Y X =∅ (#X )p X = E (# ((x 0 + X 0 ) ∩ Y)) = E   p j=1 1 (y j ∈ x 0 + X 0 )   = p E (ν d (X)) E ν d Y ⊕ X = 0.
Let us now turn to the limit of the right-hand side of Equation ( 7) when dividing by α and letting α tend to 0. First let us show that all the terms for which k ≥ 3 will tend to 0. By induction for all k ≥ 2, the terms ∂ p-k ∂t I c φ α (0, . . . , 0) have a finite limit when α tends to 0. Besides, for all k ≥ 3,

ψ (k) α (0) = i k E (T α (a) -E(T α (a))) k ≤ E |T α (a) -E(T α (a))| k ≤ 2 k α -sk , and
σ T α = α 2 -α Var(T α (a)) ∼ α→0 Var(a) 2 α 1 2 ,
where u(α) ∼ α→0 v(α) means that u(α)/v(α) tends to 1 as α tends to 0. Using the classic notation u(α) = O α→0 v(α) (meaning that there exists some constant Γ such that |u(α)| ≤ Γv(α) in the neighborhood of 0), one observes that for all k ≥ 3,

1 α α σ T α k ψ (k) α (0) I⊂{1,...,p} #I=k ∂ p-k ∂t I c φ α (0, . . . , 0)     X ⊂Y X =∅ i∈I 1(y i ∈ X ) i∈I c β 1(yi∈X ) p X     = O α→0 α 1 2 k-sk-1 .
But since s < 1/6, the above exponent 1 2 k -sk -1 is positive for all k ≥ 3. Hence all the terms for which k ≥ 3 tend to 0. Now for k = 2, we have ψ In addition remark that X ⊂Y X =∅ 1(y j1 ∈ X )1(y j2 ∈ X )p X = P ({y j1 , y j2 } ⊂ x 0 + X 0 ) = γ X (y j1 -y j2 ) E ν d Y ⊕ X .

Coming back to Equation [START_REF] Cao | Dead leaves model for measuring texture quality on a digital camera[END_REF], one sees that ∂ p φα ∂t1...∂tp (0, . . . , 0) admits a finite limit when α tends to 0. Noting (i) p L this finite limit, that is L = lim Hence, g α (y) -g T α (y) is the sum of two r.v. which tends to 0 in L 2 . This completes the proof.

Figure 1 :

 1 Figure 1: Addition of a transparent object. The transparency coefficient of the disc is α = 0.5.

Lemma 4 (

 4 Precision associated with the leaves layer). Let y ∈ R d and let fn (y) = α n-1 k=0 a(y, k)β k be the restriction of the sum defining f (y) to the n latest leaves which have fallen on y. Then f (y) -fn (y) ≤ Aβ n .

  (a) Original shape (b) α = 1 (occlusion) (c) α = 0.7 (d) α = 0.4 (e) α = 0.05 (f) α = 0.05, enhanced contrast

Figure 2 :

 2 Figure 2: TDL realizations with various transparency coefficients α. The RACS Xi are all obtained from the original shape 2(a) by applying a rotation of angle θ ∼ U nif (0, 2π) and a homothety of factor r ∼ U nif (0, 1), and Pa = U nif (0, 255). For α = 1, one obtains a colored dead leaves model. As soon as the leaves are transparent (α < 1), one can distinguish several layers of leaves and not only the leaves on top. For α = 0.05, the variance of the TDL model is nearly 0 (see Proposition 3). Enhancing the contrast of the image reveals the structure of the image (see 2(f)).

(a) α = 1 (

 1 b) α = 0.5 (c) α = 0.2 (d) α = 0.1 (e) α = 0.01 (f) Limit Gaussian r.f.

Figure 3 :

 3 Figure 3: From colored dead leaves to Gaussian random fields: Visual illustration of the normal convergence of the normalized TDL random fields fα -E (fα) Var (fα) α

Proposition 10 (

 10 A recurrence relation for the moments of a multivariate normal distribution). Let Y = (Y 1 , . . . , Y p ), p ≥ 2, be a normalized Gaussian vector. Then,

i 2

 2 Var(T α (a)). Besides, by induction,lim α→0 ∂ p-2 ∂t {j1,j2} c φ α (0, . . . , 0) = (i) p-2 E   p l∈{1,...,p}\{j1,j2} f G (y l )   .Hence, considering subsets I of {1, . . . , p} with two elements,

  ..,p}\{j1,j2} f G (y l )   γ X (y j1 -y j2 ) E ν d Y ⊕ X = 2 p {j1,j2}⊂{1,...,p} γ X (y j1 -y j2 ) E (ν d (X)) E   p l∈{1,...,p}\{j1,j2} f G (y l )   .This is exactly the recursive formula for the moments of a Gaussian vector given by Proposition 10. Hence, proof of Lemma 11.

5. 4 Lemma 12 (.

 412 Convergence in L 2 of the difference of the normalized random fields At this stage to conclude the proof of Theorem 7 we have to demonstrate the following lemma. Convergence to 0 in L 2 of g α -g T α ). Let g α and g T α be, respectively, the normalized TDL model and the normalized TDL model with truncated colors. Then for all y ∈ R d , g α (y) -g T α (y) Since a ∈ L 2 and for all b ∈ R, |T α (b)| ≤ |b| and lim α→0 T α (b) = b, by dominated convergence lim α→0 Var(a -T α (a)) = 0, and in particular lim α→0 Var(T α (a)) = Var(a). Let y ∈ R d . Recall that Var(f α ) = α 2 -α Var(a) and Var(f T α ) = α 2 -α Var(T α (a)).One hasg α (y) -g T α (y) = f α (y) -E(a) Var(f α ) -f T α (y) -E(T α (a)) Var(f T α ) = f α (y) -E(a) Var(f α ) -f T α (y) -E(T α (a)) Var(f α ) + f T α (y) -E(T α (a)) Var(f α ) -f T α (y) -E(T α (a)) Var(f T α ) = f α (y) -f T α (y) -E(a -T α (a)) Var(f α )Let us denote by I 1 (α) and I 2 (α) the two terms above. Remark that the numerator of I 1 (α) is a TDL model with color distribution a -T α (a) -E(a -T α (a)). Hence we haveE(I 1 (α) 2 ) =α 2-α Var(a -T α (a)) α 2-α Var(a) = Var(a -T α (a)) Var(a)

Acknowledgments

This work was supported by ANR MATAIM. The authors are grateful to Professors Jean-Michel Morel and François Roueff for their insightful suggestions.


Marie-Pierre Béal
email: beal@litp.ibp.fr

Olivier Carton
email: carton@univ-mlv.fr

Christophe Reutenauer

Marie-Pierre B Eal

Cyclic languages and Strongly cyclic languages

Introduction

Cyclic languages and strongly cyclic languages are two classes of languages of nite words over a nite alphabet. A cyclic language is conjugation-closed and for any two words having a power in common, if one of them is in the language, then so is the other. A strongly cyclic language is the set of words stabilizing a subset of the set of states of a nite deterministic automaton, the subset stabilized depending on the word stabilizing it. One says that the language stabilizes the automaton. A strongly cyclic language is rational.

We prove that a rational cyclic language is a boolean combination of strongly cyclic languages. More precisely, each rational cyclic language can be written as a chain of strongly cyclic languages. This result allows us to extend the computation of the zeta function (and generalized zeta function) of strongly cyclic languages done in B ea95] to the class of rational cyclic languages. The zeta function of a formal language L is (L) = exp( P a n z n n ), where a n is the number of words of length n in L. The motivations of this de nition and the connections with algebraic geometry and symbolic dynamics are discussed in BR90]. The zeta function of a strongly cyclic language L is equal to the zeta function of the so c system de ned by the nite automaton stabilized by L. The rationality and computability of the zeta function of a so c system have been established in Bow78] and Man71]. The formula of computation given in Bow78] are proved in B ea95] by the use of a construction on nite automata called external power. The rationality of the zeta function of a rational cyclic language has been established in BR90]. The result we give here leads to another proof and to a di erent computation.

We assume that the reader knows the basics of formal languages (see Eil72]). We also assume that the reader is familiar with the elementary notions of semigroup theory. For example, notions like syntactic monoid, Green relations, regular D-classes, minimal ideal and 0-minimal ideal are supposed to be known. We refer to Lal79] and Pin86] for a presentation of this subject.

The paper is organized as follows. Section 2 and 3 give the basic properties of cyclic languages and strongly cyclic languages. The chain-decomposition of a rational cyclic language in strongly cyclic languages is established in section 4. The computation of the zeta function and the generalized zeta function is done in the last section.

Cyclic languages

In this section, we introduce cyclic languages and give some basic properties. In the following, we denote by A a nite alphabet. In the sequel, M will always denote a nite monoid. Every element s of M has a power which is an idempotent. We denote by s ! this idempotent.

De nition 1. A language L of A is said to be cyclic if it satis es 8u 2 A ; 8n > 0 u 2 L , u n 2 L 8u; v 2 A uv 2 L , vu 2 L A language is cyclic if it is closed under conjugation, power and root. If L is a submonoid of A which is cyclic, it is then pure BP84].

Example 1. If A = fa; bg, the language L = A aA = A b of words having at least one a is cyclic.

Example 2. The language L = fa p b n1 a n2 b n2 : : : a n k b n k a q j n i 0 and p + q = n 1 g fb p a n1 b n1 a n2 b n2 : : : a n k b q j n i 0 and p + q = n k g is cyclic but not rational.

We will now only consider rational cyclic languages. Rational cyclic languages have the following straightforward characterization in terms of nite monoids.

Proposition 2. Let L A be a rational language. Let ' : A ! !M be a morphism from A onto a monoid M such that L = ' 1 (P ). The language L is cyclic if and only if 8s 2 M; 8n > 0 s 2 P , s n 2 P 8s; t 2 M st 2 P , ts 2 P From the previous characterization, we deduce some useful facts about the structure of the image of a rational cyclic language onto a nite monoid recognizing this language. We also deduce a property of the syntactic monoid of a cyclic language.

Corollary 3. Let ' : A ! !M be a morphism from A onto a monoid M such that L = ' 1 (P ). Let H a regular H-class of M. One has H P or H \P = ;. Proof. Let us suppose h 1 2 H belongs to P. For any h 2 2 H, we have h ! 1 = h ! 2 = e where e is the idempotent of H. We then have h ! 1 = h ! 2 2 P and h 2 2 P. Corollary 4. Let ' : A ! !M be a morphism from A onto a monoid M such that L = ' 1 (P ). Let H 1 and H 2 be two regular H-classes of a regular D-class.

If one has H 1 P, one also has H 2 P.

Proof. Let e 1 and e 2 be the respective idempotents of H 1 and H 2 . As two idempotents of a same D-class are conjugated (see Pin86]), there are two elements x 1 and x 2 of M such that x 1 x 2 = e 1 and x 2 x 1 = e 2 . If H 1 P, we have e 1 = x 1 x 2 2 P and then e 2 = x 2 x 1 2 P and the H-class H 2 satis es H 2 P by the previous corollary.

Corollary 5. The syntactic monoid of a cyclic language has a zero.

Proof. Let M be the syntactic monoid of a rational cyclic language and let D the minimal ideal of M. Then D is a completely regular D-class. By the previous corollary, we have D P or D \ P = ;. Let s be an element of D. For any x; y 2 M, we have xsy 2 D because D is the minimal ideal of M. If we have D P, the contexts of s are M M. If we have D \ P = ;, the contexts of s are ; ;. In both cases, all the elements of D are equivalent for the Nerode congruence, and D has only one element. The syntactic monoid of L has then a zero.

Strongly cyclic languages

We now de ne the notion of a strongly cyclic language.

De nition 6. Let A = (Q; A; E) be a deterministic nite automaton where Q is the set of states and E the set of transitions. We say that a word w stabilizes a subset P Q of states if we have P:w = P. This means 8p 2 P p:w 2 P 8p 0 2 P 9p 2 P p:w = p 0 We denote by Stab(A) the set of the words w such that w stabilizes a subset P of states in the automaton A. It should be noticed that in this de nition the subset P of states stabilized by w may depend on w. We point out that the empty word " stabilizes the set Q and therefore belongs to Stab(A) for any automaton A. We say that a language L is strongly cyclic if there is an automaton A such that L = Stab(A). In this case, we say that the language L stabilizes the automaton A. The terminology is justi ed since strongly cyclic languages are cyclic. It could be proved directly but we shall obtain this fact as a consequence of the characterization of strongly cyclic languages. The following result gives a characterization of the words w stabilizing a subset of states in an automaton.

Proposition 7. Let A = (Q; A; E) be a deterministic nite automaton. A word w stabilizes a subset P of states in A if and only if there is some state q of A such that for any integer n, the transition q:w n exists.

Proof. Suppose rst that the word w stabilizes the subset P of states. By denition, for any state p 2 P, the transition p:w n exists.

Conversely, suppose all the transitions q:w n exist for some state q of A. Since the automaton is nite, there are two integers l < m such that q:w l = q:w m . Let P be the set fq:w i j l i mg. It is straightforward that the word w stabilizes the subset P.

The following theorem gives a characterization of the strongly cyclic languages.

Theorem 8. Let L be a rational language di erent from A . The following conditions are equivalent.

1. The language L is strongly cyclic. 2. There is a morphism ' from A onto a monoid M having a zero such that L = ' 1 (fs 2 M j s ! 6 = 0g).

3. The syntactic monoid M(L) of L has a zero and the image of L in M(L) is fs 2 M j s ! 6 = 0g. Proof. Suppose rst that the language stabilizes the automaton A. Let M the transition monoid of A and ' the canonical morphism from A onto M. We show then that this monoid has a zero and that the image of L is equal to fs 2 M j s ! 6 = 0g. Let w be a word not belonging to L. By Proposition 7, for each q 2 Q, there is an integer n q such that the transition q:w nq does not exist.

For n greater than every n q , the transitions q:w n do not exist for any q. The transition induced by w n is then the empty transition and the element '(w n ) is a zero of the monoid M. This means that '(w) ! = 0. On the contrary, for any word w 2 L, the transition induced by w n is not the empty transition since there is a state q 2 Q such that the transition q:w n exists. This proves that '(w) ! 6 = 0. We have then proved that L is equal to ' 1 (fs 2 M j s ! 6 = 0g) Suppose now there is a morphism ' from A onto a monoid N having a zero such that L = ' 1 (ft 2 N j t ! 6 = 0g). Since the morphism is onto, the syntactic monoid M of L is a quotient of N: there is morphism : N! !M from N onto M.

The image (0) of the zero of N is then a zero of M. Since the zero of N does not belong to the image of L, the zero of M does not belong to the image of L. Let t a element of N such that t ! = 0. We have (t) ! = (t ! ) = 0. On the contrary, if t ! 6 = 0, the element t ! belongs to the image of L, and so does (t ! ). This implies (t ! ) 6 = 0. Finally, we have ft 2 N j t ! 6 = 0g = 1 (fs 2 M j s ! 6 = 0g).

Suppose that the syntactic monoid M of L has a zero and that the image of L in M is fs 2 M j s ! 6 = 0g. We denote by ' the canonical morphism from A onto M. We build the following deterministic automaton A = (Q; A; E). The set of states of A is the set Q = M f0g of the non-zero elements of M. The transition q:a is q:a = q'(a) if q'(a) 6 = 0 and does not exist otherwise. It can be easily checked that a transition q:w is q:w = q'(w) if q'(w) 6 = 0 and does not exist otherwise. We show now that the language L stabilizes the automaton A.

Let w a word in L. For q = '(w), the transition q:w n is q:w n = '(w) n+1 and exists since '(w) n+1 6 = 0. On the contrary, for a word w not in L, there is an integer m such that '(w) n = 0. The transition q:w n does not exist for any q. By the previous Proposition, the language L stabilizes the automaton L. This nishes the proof.

Corollary 9. The previous characterization shows that strongly cyclic languages are cyclic. The previous theorem can be used to prove that a given language is not strongly cyclic. Without this characterisation, such results are sometimes hard to obtain.

Example 4. The syntactic monoid of the language L = A aA is the two elements monoid M = fb = 1; a = 0g. The D-class structure of M is shown in Figure 2.

The image of L in M is the singleton f0g and the previous theorem states that the language L is not strongly cyclic.

Decomposition of cyclic languages

In this section, we prove the main result.

By the de nition of cyclic languages, a boolean combination of cyclic languages is still a cyclic language. In particular, a boolean combination of strongly cyclic languages is a rational cyclic language. The following result gives some converse.

Theorem 10. Any rational cyclic language is a boolean combination of strongly cyclic languages.

The proof of the theorem is based on the following lemma. By a strict quotient of M, we mean a quotient which is strictly smaller than M.

Lemma11. Let L be a rational cyclic language and ' : A ! !M a morphism from A onto a nite monoid M such that L = ' 1 (P ). Let suppose furthermore that the monoid M has a zero and that this zero does not belong to P.

Then, either L is recognized by a strict quotient of M or there exists a strongly cyclic language L 0 such that L L 0 and such that the language L 0 L is recognized by a strict quotient of M.

In both cases, the zero of the quotient does not belong to the image of the language (L in the rst case and L 0 L in the second case).

Proof. Let P the image of L in M. Suppose D 1 ; : : : ; D n are the 0-minimal Dclasses of M. For any s 2 D i and x; y 2 M, we have xsy = 0 or xsy 2 D i .

Suppose D 1 satis es D 1 \P = ;. For any x; y 2 M, we have xsy 6 2 P. All the elements of D 1 are equivalent to 0 by the Nerode congruence. The language L is then recognized by the Rees quotient M=I where I is the ideal I = D 1 f0g. We can suppose that every D i satis es D i \ P 6 = ;. There exists then s i 2 D i \ P. Since s i 2 P, the idempotent s ! i is not zero and belongs then to D i because this D-class is 0-minimal. The element s i belongs then to a regular Hclass H i of D i . Let s be an element of D i . If s ! = 0, the element s does not belong to P because 0 does not belong to P. On the contrary, If s ! 6 = 0, the element s belongs to a regular H-class H 0 i of D i . Since H i \ P 6 = ;, we have H i P and H 0 i P by Corollaries 3 and 4. Finally, we have proved that D i \ P = fs 2 D i j s ! 6 = 0g Let P 0 be de ned by P 0 = fs 2 M j s ! 6 = 0g. By Theorem 8, the language L 0 = ' 1 (P 0 ) is strongly cyclic. It is straightforward that P P 0 . Since we have D i \ (P 0 P) = ;, the language L 0 L is recognized by the Rees quotient M=I where the ideal is equal to f0g S n i=1 D i . This quotient is strictly smaller than M and this nishes the proof of the lemma.

The following lemma states that the class of strongly cyclic languages is closed under union and intersection.

Lemma 12. Let L 1 and L 2 be two strongly cyclic languages. Both languages L 1 L 2 and L 1 \ L 2 are then strongly cyclic.

Proof. We suppose that L 1 and L 2 respectively stabilizes the automata A 1 = (Q 1 ; A; E 1 ) and A 2 = (Q 2 ; A; E 2 ). We can suppose that Q 1 \ Q 2 = ;. The language L 1 L 2 stabilizes then the automaton A 1 A 2 . The language L 1 \ L 2 stabilizes the automaton A 1 A 2 = (Q 1 Q 2 ; A; E 3 ) where the transition E 3 is de ned by (p 1 ; p 2 ):a = (q 1 ; q 2 ) if the transitions p 1 :a = q 1 and p 2 :a = q 2 exist.

We can now complete the proof of the theorem.

Proof. We prove that every cyclic language L can be written as a chain of strongly cyclic languages. This means that there are strongly cyclic languages L 1 ; : : : ; L n satisfying L 1 L 2 L n such that

L = L 1 L 2 + L 3 L n
We prove the result by induction on the size of a nite monoid M having a zero and recognizing the language L. We suppose that there is a morphism ' :

A ! !M from A onto a nite monoid M such that L = ' 1 (P ). We also suppose that the monoid M has a zero. By Corollary 5, the syntactic monoid of L has this property. If the monoid M has only one element, the language L is either ; or A which are both strongly cyclic. The empty language ; stabilizes the empty automaton. The full language stabilizes the automaton having one state and a transition for each letter from this state to this state.

If the zero of M belongs to the image of L in M, we replace L by A L which is also cyclic. It is then su cient to prove the result for A L. Indeed, if the complement A L of L can be written as a chain A L = L 1 L 2 + L n , the language L can be written L = A L 1 + L 2 L n which is also a chain of strongly cyclic languages.

We can now suppose that the zero of M does not belong to the image P of L in M. By Lemma 11, either L is recognized by a strict quotient of M and the induction hypothesis immediately applies or there is a strongly cyclic language L 0 such that L 0 L is recognized by a strict quotient of M. By the induction hypothesis, the language L 0 L can be written as a chain of strongly cyclic languages, i.e., L 0 L = L 2 L 3 + L n . The language L can be then written as the chain L = L 1 L 2 + L n where the language L 1 is equal to L 0 L 2 which is strongly cyclic by Lemma 12.

We point out that the zero of the smaller monoid recognizing L 0 L does not belong to the image of L 0 L. It is not necessary to replace this language by its complement any more. We give here an example. Example 5. Let L be the language (b + aa) + (ab a) + a b . The structure of the syntactic monoid of L is given on gure 3. The image P of L in M(L) is equal to P = fa; aa; aba; aabg.

The subset P 0 de ned in the proof is equal to P 0 = f1; a; aa; b; aba; aabg and the language L 0 is (b + aa) + (ab a) + a . The language L 0 L is then equal to b which is strongly cyclic. The languages L 0 and L 0 L stabilize respectively the automata A 1 and A 2 (see Figure 1 p. 4).

We have directly proved that every rational cyclic language can be written as a chain of strongly cyclic languages. In fact, it is just necessary to prove that rational cyclic languages are boolean combination of strongly cyclic languages.

A general result states that if a class F of sets is closed under union and intersection, every set belonging to the boolean closure of F can be written as a chain of sets of F. For further details see Car93] (chapter 3).

Zeta function of a cyclic language

We rst give the de nitions of generalized zeta function and zeta function of a language of nite words over a nite alphabet A.

If A is a nite alphabet, we note ZhhAii (resp. Z A]]) the algebra of non commutative (resp. commutative) formal series with coe cients in Z over the alphabet A. The subset of non commutative polynomials is denoted by ZhAi. We note ' the natural algebra homomorphism from ZhhAii to Z A]] which makes the letters commute. For example, '(2ab 3ba) = ab. We note the morphism from Z A]] to Z z]] de ned by (a) = z for each letter a of A.

Fig. 1 .

 1 Fig. 1. Automata A1 and A2

Fig. 2 .

 2 Fig. 2. Structure of the syntactic monoid of L.

Fig. 3 .

 3 Fig. 3. Structure of the syntactic monoid of L.

Let L be a language of nite words over A, we note L the characteristic series of L. This series belongs to ZhhAii and admits the decomposition:

where L n is the homogeneous part of degree n of L.

De nition 13. The generalized zeta function of a language L over the alphabet A is now the following commutative series:

De nition 14. The zeta function of a language L over the alphabet A is the following series :

where a n is the number of words of L of length n.

It is shown for example in B ea95] that the generalized zeta function of a strongly cyclic language L of an automaton A is the generalized zeta function of the so c system de ned by A, that is the set of bi-in nite words that are labels of bi-in nite paths of A. The zeta function of a so c system counts periodic orbits of the symbolic dynamic system. The zeta function of a strongly cyclic language is a rational series. The computation of the generalized zeta function and the zeta function of a strongly cyclic language was done in B ea95] by using a construction on nite automata called external power.

If A = (Q = f1; 2; : : :; ng; E; T) is a deterministic automaton, the external power of order k, where 1 k jQj, of the automaton A is the automaton (Q 0 ; E 0 ) labelled in f 1; 1g A, where Q 0 = f(i 1 ; i 2 ; :::; i k ) 1 i1<i2< <i k n g.

There exists an edge labelled ( )a from (i 1 ; i 2 ; :::; i k ) to (j 1 ; j 2 ; :::; j k ) if for each l with 1 l k, there exists one edge in E labelled a going out from i l and (j 1 ; j 2 ; :::; j k ) is the image of (i 1 a; : : : ; i k a) by a permutation of signature ( ). The automaton A is equal to its external power of order 1 by identi cation of +a to a. The matrix associated to an automaton labelled on f 1; 1g A is the square matrix (x ij ) 1 i;j n where x pq is the sum (in ZhAi) of the labels of edges from p to q. The commutative matrix associated is the matrix ('(x ij )) 1 i;j n .

We denote by Q i the commutative matrix associated to the external power of order i of the automaton A. We then have (see

where I is the identity matrix of the same size as Q i .

Computation of the zeta function of a cyclic language

The result of section 4 can be used to extend the previous computation of zeta functions of strongly cyclic languages to all cyclic languages. This gives an other proof of the rationality of the zeta function of a cyclic language established in BR90]. The computation is the following:

Let L be a cyclic language. By section 4 it can be written as a chain L 1 L 2 + + ( 1) r 1 L r , where L j+1 L j for 1 j (r 1) and where all L j are strongly cyclic languages. By de nition of the generalized zeta function we have:

(Z(L j )) ( 1) j 1 :

Example 6. We compute the generalized zeta function Z(L) and the zeta function (L) of the cyclic language L = (b + aa) + (ab a) + a b introduced in example 5. The language L admits the chain-decomposition L = L 1 L 2 where L 1 = Stab(A 1 ) (see gure 1) and L 2 = Stab(A 2 ) (see gure 1). We get:

Z(L 1 ) = j1 + aj 1 b a a 1 Z(L 2 ) = 1 j1 bj Z(L) = Z(L 1 ) Z(L 2 ) = (1 + a)(1 b) 1 b aa (L) = 1 z 2 1 z z 2 :


Emilios Cambouropoulos
email: emilios@mus.auth.gr

Maxime Crochemore
email: maxime.crochemore@univ-mlv.fr

Costas Iliopoulos

Manal Mohamed

Marie-France Sagot
email: .sagot@inria.fr

A PATTERN EXTRACTION ALGORITHM FOR ABSTRACT MELODIC REPRESENTATIONS THAT ALLOW PARTIAL OVERLAPPING OF INTERVALLIC CATEGORIES

Keywords: string, don't care, repetitions, suffix tree, lowest common ancestor

This paper proposes an efficient pattern extraction algorithm that can be applied on melodic sequences that are represented as strings of abstract intervallic symbols; the melodic representation introduces special "don't care" symbols for intervals that may belong to two partially overlapping intervallic categories. As a special case the well established "step-leap" representation is examined. In the step-leap representation, each melodic diatonic interval is classified as a step (±s), a leap (±l) or a unison (u). Binary don't care symbols are introduced to represent the possible overlapping between the various abstract categories e.g. * = s, * = l and # = -s, # = -l. For such a sequence, we are interested in finding maximal repeating pairs and repetitions with a hole (two matching subsequences separated with an intervening non-matching symbol). We propose an O(n + d(nd) + z)-time algorithm for computing all such repetitions in a given sequence x = x[1..n] with d binary don't care symbols, where z is the output size.

INTRODUCTION

Recently, there have been different proposals in the literature to develop an effective music information retrieval system. The goal of these proposals is to take advantage of appropriate computer science techniques. For example, representing the musical surface as a string or set of strings may make it possible in some cases to apply existing algorithms from the field of stringology. For instance, in order to discover similarities between different musical entities or to establish motivic "signatures", music analysts may use algorithms that extract repetitions from strings. Such similarities often involve finding approx-Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. c 2005 Queen Mary, University of London imate repetitions [START_REF] Crawford | String matching techniques for musical similarity and melodic recognition[END_REF]. This requires developing new approximation measures that meet musicians' needs.

One commonly used representation for music is the numeric representation MIDI. For such a representation, different approximation measures have been developed, such as, δ-, γand {δ, γ}-approximate. For example, in δ-approximate matching, equal-length strings consisting of integers match if each corresponding integer differs by not more than δ -e.g. a C-major {60, 64, 65, 67} and a Cminor {60, 63, 65, 67} sequence can be matched if a tolerance δ = 1 is allowed in the matching process. Using these approximation measures, algorithms for finding approximate repetitions in musical sequences have been developed [START_REF] Iliopoulos | Computing approximate repetitions in musical sequences[END_REF][START_REF] Cambouropoulos | Algorithms for computing approximate repetitions in musical sequences[END_REF]. These algorithms are based on approximate pattern matching techniques. For an overview refer to Clifford and [START_REF] Clifford | Approximate string matching for music analusis[END_REF].

Although MIDI is the most common representation in the computational domain, it has certain well-known shortcomings, for instance, many important musical properties are not explicitly represented (e.g. note durations, accidentals etc.) and almost all information on musical structure is lost. Therefore, different representations have been proposed in the literature. For example, [START_REF] Hawley | Structure of Sound[END_REF] proposed representing the musical signal as a sequence of pitch intervals. In order to allow tolerance in interval matching, [START_REF] Ghias | Query by humming: Musical information retrieval in an audio database[END_REF] used the reduced interval alphabet of the "melodic contour" representation. [START_REF] Lemström | Musical information retrieval using musical parameters[END_REF] proposed classifying the intervals into seven partially overlapping classes: small, medium and large, up-or downwards, and prime.

In this paper, we propose an alternative method to using approximate pattern matching techniques for finding approximate repetitions in a musical string. Our approach is based on using exact pattern matching techniques to extract repetitions from an abstract level of a musical sequence. As an abstract representation, we will use the "refined contour" (or step-leap) representation -see, for instance, application of this representation in <http://www.themefinder.org>. In the step-leap representation, intervals are classified into five distinct equivalence classes: up-or downwards step and leap, and unison. An interval with magnitude a = 0 is a unison (u), a < 2 is a step (s), and any other interval a ≥ 2 is a leap (l); the direction of intervals is preserved -see second In the second string of symbols in Figure 1, two repeated substrings are found: -s s -l and l -s s s -l each occurring twice. However, for a listener/musician, the second half of this string of intervals is an approximate repetition of the first half (two approximately matching substrings separated by a "hole" of size one) because intervals a = 1 and a = 2 are considered similar (i.e. a step is similar to a small leap). This is not simply some rare exception in music. It is a rather common phenomenon especially when themes appear in their dominant form (see, for instance, the tonal answers of almost half of Bach's fugue themes from the two books of the Well-Tempered Clavier). In Figures 2 and3 some melodic examples are presented.

The problem in the step-leap representation is that the abstract interval classes (u, s, l) have sharp boundaries and no diatonic pitch interval instance may belong to more than one class. In other words, borderline members can never be matched to other 'similar' members of other classes (e.g. an a = 2 interval as a member of leap can never be matched to a 'similar' a = 1 interval which is a step), i.e, a small leap can never be considered as a step. A way to overcome this problem is to allow partial overlapping between the various classes (this is also suggested in [START_REF] Lemström | Musical information retrieval using musical parameters[END_REF]). For instance, an interval a = 2 may be classified as either step or leap. A special "binary don't care" symbol * that matches either s or l is used (see third string of symbols in Figure 1). Similarly, another don't care symbol # that matches both -s and -l is also used. Note that this idea can be easily extended to any constant number of partially overlapping classes, such as those proposed in [START_REF] Lemström | Musical information retrieval using musical parameters[END_REF].

A pattern extraction algorithm may find a large number of repeating patterns. Many of these are not musically or perceptually important. A mechanism, therefore, for selecting important patterns is required. A relatively sophisticated method for finding beginnings of potentially salient repetitions is proposed by Cambouropoulos. Following this proposal, for each pattern a prominence value is calculated based on frequency of occurrence, pattern length and degree of pattern overlapping; these prominence values contribute to establishing a segmentation prominence profile for a melody whereby the most likely positions of important repetitions are highlighted. In this paper, where the aim is not segmentation but the extraction of interesting patterns, simpler criteria are set: extraction of maximal repeating pairs and repetitions with a 'hole'. The later involves extracting immediately repeated nonoverlapping melodic patterns -in order that the patterns do not overlap over one note a 'hole' is necessary in the pitch interval representation (for example, in Figure 1, the second half of the melody is an approximate repetition of the first half -in the interval representations beneath the melody it is necessary to skip one symbol in the middle so that the repetition is consecutive and non-overlapping at the note level).

For strings with don't cares, several string matching algorithms have been proposed (see [START_REF] Fischer | String matching and other products[END_REF][START_REF] Apostolico | Optimal off-line detection of repetitions in a string[END_REF][START_REF] Amir | Approximate subset matching with don't cares[END_REF]). Recently, [START_REF] Iliopoulos | String regularities with don't cares[END_REF] presented algorithms for computing typical regularities in strings with don't cares. Here, we consider binary don't care symbols that each matches beside itself two additional symbols. For string x = x[1..n] with d binary don't cares, we propose an algorithm for computing special kinds of repetition that we refer to as "maximal-pairs" and "repetitions with a hole". The proposed algorithm uses O(n + d(nd) + z) time, where z is the output size.

The paper is organized as follows: in Section 2, we state the preliminaries used throughout the paper. In Section 3, we define all approximate repetitions problem and describe in general how to find them. In Section 4, we detail our algorithm. Finally, in Section 5, we analyze the running time of the algorithm.

PRELIMINARIES

Throughout the paper, x = x[1..n] denotes a string of length n over Σ ∪ { * , #}, where Σ = {s, -s, l, -l, u}. The symbols ' * ' and '#' are called binary don't care symbols. Each binary don't care symbol matches itself and two different symbols, that is, * = * , * = s, * = l, # = #, # = -s and # = -l.

We use x[i], for i = 1, 2, . . . , n, to denote the i-th symbol of x, and x[i..j] as a notation for the substring

x[i]x[i + 1] • • • x[j] of x. If x = uv then x is said to be the concatenation of the two strings u and v. A string y is said to occur in x at position i if y[j] = x[i + j -1], for 1 ≤ j ≤ |y|.
A repeating pair in x is represented by (p; i, j) where, x[i..i+p-1] = x[j..j+p-1] for some i = j. The positive integer p is called the period of the repeating pair. If

x[i - 1] = x[j -1] then (p; i, j) is left-maximal.
[i + p] = x[j + p] then (p; i, j) is right-maximal. If (p; i, j
) is both left-and right-maximal then it is maximalpair. A repetition with a hole is a repeating pair (p; i, j) such that j = i + p + 1.

Here, we present a method for finding all maximalpairs and all repetitions with a hole in a given string x, where x may have occurrences of binary don't cares. Our method uses the suffix tree of x as a fundamental data structure. A complete description of suffix trees is beyond the scope of this paper, and can be found in [START_REF] Gusfield | Algorithms on Strings, Trees and Sequences: Computer Science and Computational Biology[END_REF] or [START_REF] Crochemore | Jewels of Stringology[END_REF]. However, for the sake of completeness, we will briefly review the notion. Definition 1 (Suffix tree) A suffix tree T (x) of the string x$ = x[1..n]$ is a rooted directed tree with exactly n leaves numbered 1 to n + 1, where $ / ∈ Σ. Each internal node, other than the root, has at least two children and each edge is labelled with a non-empty substring of x. No two edges out a node can have edge-labels beginning with the same symbol. The key feature of the suffix tree is that for any leaf i, the concatenation of the edge-labels on the path from the root to leaf i exactly spells out i-th suffix of x, with n + 1 denotes the empty suffix.

Several algorithms construct the suffix tree T (x) in Θ(n) time and space, assuming constant size alphabet (see for example [START_REF] Crochemore | Jewels of Stringology[END_REF] and (Gusfield, 1997)). For any node v, the path-label of v is the label of the path from the root of T (x) to v; it is denoted by label(v). The string-depth of v is the number of symbols in v's path-label; it is denoted by depth(v). The leaflist of v is the set of the leaf numbers in the subtree rooted at v; it is denoted by LL(v).

Our method makes use of the Schieber and Vishkin (1988)'s Lowest Common Ancestor algorithm. For a given rooted tree T , the lowest common ancestor (LCA) of two node u and v is the deepest node in T that is ancestor of both u and v. After a linear amount of preprocessing of a rooted tree, any two nodes can be specified and their lowest common ancestor is found in constant time. That is, a rooted tree with n nodes is first preprocessed in O(n) time, and thereafter any lowest common ancestor query takes only a constant time to be solved, independent of n.

In the context of suffix trees, the situation commonly arises that both u and v are leaves in T (x), where x[i..n] and x[j..n] are the suffixes represented by u and v respectively, for integers i and j in the range 1..n + 1. In this case, the node w = LCA(u, v) is the root of the minimum size subtree contains u and v. Note that the path-label of w (label(w)) is the longest common prefix of x[i..n] and x[j..n]. The ability of finding such longest common prefix is an important primitive in many string problems.

FINDING ALL REPETITIONS PROBLEM

Here, we study the problem of finding in a given sting x over Σ ∪ { * , #}, the following two kinds of repetitions: all maximal-pairs and all repetitions with a hole, where each repetition with a hole is a repeating pair (not necessary maximal) in which an intervening symbol separates the two matching substrings. In fact, there is a very close relation between these two kinds of repetitions. For example, if a given string x is as follows: x =s s l s s # l s * -l l s s

Then the maximal-pair (6; 3, 7) -which represents two overlapping matching substrings-represents four repetitions with a hole: (3;3,7),(3;4,8), (3;5,9), (3;6,10).

The relation between maximal-pairs and repetitions with a hole is presented in the following lemma:

Lemma 1 Let (p; i, j) be a maximal-pair in x and let g = jip.

1. If g = 1 then (p; i, j) is a repetition with a hole.

2. If g ≤ 0 and j = i + 1 then (p ′ ; i + k, j + k) are repetitions with a hole, where p ′ = ji -1 and k = 0..|g| + 1.

3. If g ≤ 0 and j = i + 1 then some repetitions with a hole may be found in x[i..j + p -1].

4. If g > 1 then x[i..i + p -1] and x[j..j + p -1]g are two matching substrings separated by intervening g symbols .

In the following we will introduce Gusfield's algorithm for finding all maximal-pairs in a given string without don't cares. The basic tool behind Gusfield's algorithm is the suffix tree. The algorithm starts by constructing the suffix tree for a given string. The algorithm then uses a bottom-up approach (from leaves to root) to report for each internal node the maximal-pairs associated with it. This is accomplished by maintaining the leaf list LL(v) of each internal node v as a collection of disjoint sublists LL α (v), where α is the symbol preceding the suffix associated to a leaf in the subtree rooted by v. Thus, each internal node is attached at most |Σ| sublists. Reporting the maximal-pairs is accomplished by the cartesian product of a leaf-sublist with all the leaf-sublists of its brothers that correspond to different symbols. The algorithm runs in O(n + z), where z is the number of reported maximalpairs. Gusfield's algorithm can be easily modified to find all repeating pairs (not necessary maximal) in x.

The presence of the binary don't care symbols in x complicates the construction of the suffix tree. Hence, Gusfield's algorithm cannot be used directly to solve the problem. The dynamic programming seems to be an obvious solution. The cost of this method is quadratic. For example, if x = s s # l l * -l s * -l l l s. Then using dynamic programming, the following maximal-pairs can be found: (1;1,2), (1;1,6), (6;1,8), (1;1,9), (1;1,13), (2;2,6), (1;2,8), (1;2,13), (2;4,5), (1;4,6), ..., (1;11,12). Note that only (6;1,8), (1;4,6), (1;6,8), (1;9,11) and (1;4,6) are the only repetitions with a hole in x (see Table 1).

In the next section, we will explain how the suffix tree can be used to speed up the dynamic programming calculations. Independently of the size of the alphabet, our algorithm works for any string that have occurrences of finite number of binary don't cares.

s # l l * -l s * -l l l s 1 s - 1 0 0 0 1 0 1 1 0 0 0 1 2 s - 0 0 0 1 0 1 2 0 0 0 1 3 # - 0 0 0 2 0 0 3 0 0 0 4 l - 1 1 0 0 1 0 4 1 0 5 l - 2 0 0 1 0 1 5 0 6 * - 0 1 1 0 1 2 6 7 -l - 0 0 2 0 0 0 8 s - 1 0 0 0 1 9 * - 0 1 1 1 10 -l - 0 0 0 11 l - 1 0 12 l - 0 13 s -

ALGORITHM

Given a string x over Σ ∪ { * , #}, we construct two new strings x s and x l . Where string x s (respectively, x l ) is obtained by substituting each * by s and # by -s (respectively, each * by l and # by -l). The idea is to construct two strings both over Σ each is a complement of the other in a sense that for each binary don't care symbol in the original string each of the two new constructed strings contains one of the two possible matching symbols. Note that the suffix trees of the two constructed strings can be built in linear time. Given x s and x l , each maximal-pair (p; i, j) in x can be considered as the concatenations of m right-maximal repeating pairs:

(p 1 ; i, j), (p 2 ; i+p 1 , j+p 1 ), ..., (p m ; i+

m-1 k=0 p k , j+ m-1 k=1 p k ),
where 1. the starting-pair (p 1 ; i, j) is a maximal-pair (i.e. leftand right-maximal) in either x s or x l , 2. the collection of these right-maximal repeating pairs is distributed between x s and x l i.e. one rightmaximal repeating pair is in x s and the following repeating pair is in x l , 3. p = m k=1 p k . The above states the main idea of our algorithm. The algorithm iterates twice. In the first iteration, all maximalpairs in x whose starting-pairs are in x s are calculated. In the second iteration, all maximal-pairs in x whose starting-pairs are in x l are calculated.

Recall that the starting-pair needs to be maximal. Thus, each iteration starts by calculating all maximalpairs using the suffix tree (as in Gusfield). Then, each maximal-pair is extended to the right by a sequence of right-maximal pairs using a series of jumps from one suffix tree to another. In each attempt of jump we calculate the depth of the lowest common ancestor of two nodes. For example, if the starting-pair (p 1 ; i, j) is a maximalpair in x s then p 2 is equal to the depth of lowest common ancestor of the two leaves i + p 1 and j + p 1 in T (x l ).

Similarly, p 3 is the depth of the lowest common ancestor of leaves i + p 1 + p 2 and j + p 1 + p 2 in T (x s ) and so on.

For example, if x = s s # l l * -l s * -l l l s then x s = s ss l l sl s sl l l s and x l = s sl l l ll s ll l l s. The suffix trees of x s and x l are represented in Figures 4 and5.

Consider node u 1 ∈ T (x s ). During the bottom-up traversal of T (x s ) and at node v 1 , the maximal-pair (2;1,8) is calculated. To check whether this starting-pair can be extended to the right, that is whether x[1 + 2] matches x[8 + 2]. Since they match, the algorithm jumps to T (x l ) and calculates the lowest common ancestor of leaves 1+2 and 8+2. The lowest common ancestor of these two leaves is v 2 . Since depth(v 2 ) = 3, the current repeating pair is extended to the right by the right-maximal repeating pair (3;3,10). Since x[3 + 3] matches x[10 + 3], the algorithm jumps back to T (x s ) calculating the lowest common ancestor of the two leaves 3+3 and 10 +3, that is v 3 . Since depth(v 3 ) = 1, the current repeating pair is extended further to the right by the right-maximal repeating pair (1;6,13). Because x[6 + 1] does not match x[13 + 1], no more jumps are possible. Then, the algorithm reports (6;1,8) as a maximal-pair in x. Moreover, since g = 8 -1 -6 = 1, then (6;1,8) is a repetition with a hole in x (Lemma 1).

The details of the algorithm are presented in Figures 6 and7. For simplicity, algorithm Find-Maximal-Pairs assumes that both T (x s ) and T (x l ) are binary suffix trees. This is always a valid assumption since that any suffix tree can be transformed into binary one in O(n) time. The function Report(p; i, j) reports the given maximalpair and checks according to Lemma 1 for all repetitions with a hole. Note that if g = jip ≤ 0 and j = i + 1 then All-Pairs is used to find all pairs (not necessary maximal) in x[i..j + p -1] . Algorithm All-Pairs is a simple modification of All-Maximal-Pairs.

RUNNING TIME

In this section, we analyze the running time of All-Repetitions algorithm. Recall that, for constant size alpha- bet, a suffix tree can be built in O(n)-time. Thus, creating both T (x s ) and T (x l ) costs O(n)-time. Creating the leaflists of all leaves costs O(n)-time. At every internal node, the algorithm reports the repetitions associated with this node and constructs the leaf-sublists by concatenating the leaf-sublists of the children of this node. The total cost for creating the leaf-lists over all internal node is O(n)-time.

For each possible starting-pair the algorithm performs series of jumps form one tree to another. Each jump costs constant time which is the cost of the lowest common ancestor (LCA) query [START_REF] Schieber | On finding lowest common ancestors: Simplification and parallelization[END_REF]. In the following we will estimate an upper bound for the number of jumps performed by the algorithm.

Observe that, we jump from T (x s ) to T (x l ) to extend the current repeating pair to the right by a right-maximal repeating pair in x l . This only possible if and only if the first two symbols of both two copies of this new rightmaximal repeating pair are either * and l or # and -l. Similarly, we jump back from T (x l ) to T (x s ) if and only if the current repeating pair can be extended to the right by

. for i = 1 to n 2. if x[i] = ' * ' 3. then x s [i] = 's' 4. else if x[i] = '#' 5. then x s [i] = '-s' 6. else x s [i] = x[i] 7. for i = 1 to n 8. if x[i] = ' * ' 9. then x l [i] = 'l' 10. else if x[i] = '#' 11. then x l [i] = '-l' 12. else x l [i] = x[i] 13.
Build the suffix trees T (x s ) and T (x l ) 14. Find-Maximal-Pairs(x, T (x s )) 15. Find-Maximal-Pairs(x, T (x l ))

Figure 6: All-Repetitions algorithm Algorithm Jump&Report(x, i, j, c, d) 1. length ←d 2. while x[i + length] = x[j + length] 3. if c = 's' then c ←'l' 4.
else c ←'s' 5.

v ←T (x c ).LCA(i + length, j + length) 6.

length ←length + depth(v) 7. Report(length; i, j)

Algorithm Find-Maximal-Pairs(x, T (x c )) 1. for each leaf node u ∈ T (x c ) 2. if u represents the ith suffix of x c 3. then LL x[i-1] (u) ←{i} 4. for each α ∈ Σ ∪ { * , #} and α = x[i -1] 5.
LL α (u) ←∅ 6. for each internal node u ∈ (x c ) in bottom-up (depthfirst) manner 7.

u 1 , u 2 ←the left and the right children of u 8.

for each (i ∈ LL α1 (u 1 ) and j ∈ LL α2 (u 2 ))

where

α 1 = α 2 9. if (x[i + depth(u)] = x[j + depth(u)]) 10.
then Jump&Report(x,i,j,c,depth(u)) 11.

else Report(depth(u); i, j) 12.

for each α ∈ Σ ∪ { * , #} 13.

LL α (u) ←LL α (u 1 ) ∪ LL α (u 2 ) Clearly, the algorithm might have a quadratic running time if the input string has n/2 binary don't care symbols. For example, finding all repetitions in string x = {sl} n/4 * n/2 will cost O(n 2 )-time. This is asymptotically equal to the running time of the dynamic programming. In practice, we expect our algorithm to have a better performance. Table 2 shows the values in the dynamic programming matrix that are calculated using All-Repetitions algorithm to compute all maximal-pairs and all repetitions with a hole in string x = s s # l l * -l s * -l l l s. Note that, in addition to the 22 reported repetitions, only 4 intermediate values have been calculated by our algorithm. In this paper we have presented an algorithm that enables extraction of melodic patterns from abstract strings of symbols; this abstract representation allows partial overlapping between the various abstract symbolic classes. As a special case, we have applied the proposed algorithm on the commonly used "step-leap" interval representation.

In terms of melodic representation, it is suggested that a more refined representation that comprises of a larger number of abstract interval classes (e.g. unison, step, small leap, medium leap, large leap) may actually enable the extraction of better melodic patterns from the standpoint of musical analysis or, even, music perception. Additionally, the use of rhythmic 'contour', in terms of rhythmic abstract classes (e.g. equal, slightly larger, larger, much larger), may improve results further. Such representations have yet to be studied, implemented and tested.

The proposed algorithm requires extensive testing on pattern extraction tasks, and its performance has yet to be compared with other similar algorithms. This study, however, has presented a novel problem in terms of melodic representation and pattern extraction, and has attempted to provide an efficient solution to it that can be used for further testing and evaluation.

Figure 1 :

 1 Figure 1: Melodic pattern-matching example (the pitches of this example are taken from Bach's Well-Tempered Clavier, Book I, Fugue in F# major)

  Figure 2: Melodic examples where an a = 1 interval (step) and an a = 2 (leap) should be matched (these positions are indicated by asterisks in the melodic examples). Brackets indicate extracted melodic repetitions

Figure 5 :

 5 Figure 4: The suffix tree of x s $

Table 1 :

 1 Using dynamic programming to find all repetitions. The bold values represent the lengths of all maximal-pairs

	1	2	3	4	5	6	7	8	9	10 11 12 13
	s

  Figure7: Jump&Report and Find-Maximal-Pairs subroutines a right-maximal repeating pair in x s , where the first two symbols of both copies of this repeating pair are either * and s or # and -s. Thus, the total number of jumps is O(d(nd)), where d is the total number of don't cares in x. Summing the above gives that the total running time is as follows: Given string x[1..n] ∈ {Σ ∪ { * , #}} * , algorithm All-Repetitions reports all maximal-pairs and all repetitions with a hole in x in space O(n) and timeO(n + z + d(nd)), where z is the output size and d is the total number of binary don't cares.

	Theorem 1

Table 2 :

 2 The values calculated and reported by All-Repetitions algorithm. The bold values represent the lengths of the reported maximal-pairs

s -6 CONCLUSIONS

ACKNOWLEDGEMENTS

Maxime Crochemore is partially supported by CNRS, Wellcome Foundation, and Nato grants. Costas S. Iliopoulos is partially supported by a Marie Curie fellowship, Wellcome Foundation, Nato and Royal Society grants. Manal Mohamed is supported by an EPSRC studentship. Marie-France Sagot is partially supported by French Programme BioInformatique Inter EPST, Wellcome Foundation, Royal Society and Nato grants.