Toni Hotanen
email: tonhot@utu.fi

Everywhere Zero Pointwise Lyapunov Exponents for Sensitive Cellular Automata

Keywords: cellular automata, sensitive, Lyapunov exponents

Lyapunov exponents are an important concept in differentiable dynamical systems and they measure stability or sensitivity in the system. Their analogues for cellular automata were proposed by Shereshevsky and since then they have been further developed and studied.

In this paper we focus on a conjecture claiming that there does not exist such a sensitive cellular automaton, that would have both the right and the left pointwise Lyapunov exponents taking the value zero, for each configuration. In this paper we prove this conjecture false by constructing such a cellular automaton, using aperiodic, complete Turing machines as a building block.

Introduction

Cellular automata are discrete dynamical systems, consisting of a set of configurations over a regular lattice of cells, a finite set of symbols, a neighbourhood vector of cells and a local update rule. The local update rule, together with the neighbourhood, defines a global rule. The global rule takes each cell to a state dictated by the local rule, when given the states of a given cell's neighbours as an input. The regular lattice we are interested in, is often a countable group. Classically the group is chosen to be Z d for some d ∈ Z + , where d is called the dimension of the cellular automaton. The set of configurations is the set of all mappings from the selected group to the finite set of symbols. The global function is continuous with respect to the prodiscrete topology equipped to the set of configurations.

Lyapunov exponents are a measure for the rate of divergence of infinitesimally close trajectories in dynamical systems. They were first introduced in 1892 by Lyapunov in his doctoral thesis titled: The general problem of the stability of motion, English translation of which can be found in [START_REF] Lyapunov | General Problem of the Stability Of Motion[END_REF]. Since then they have been widely studied in the context of differentiable dynamical systems. Their importance in the study of non-linear dynamical systems, for example, can be found stated in [START_REF] Greiner | Lyapunov Exponents and Chaos[END_REF]. In the context of cellular automata, the Lyapunov exponents were first considered by Wolfram in [START_REF] Wolfram | Universality and complexity in cellular automata[END_REF]. In Problem 2 of [START_REF] Wolfram | Twenty problems in the theory of cellular automata[END_REF], the question to establish the exact connection between entropies and Lyapunov exponents is asked. A first formal definition of Lyapunov exponents for one-dimensional cellular automata is due to Shereshevsky in [START_REF] Shereshevsky | Lyapunov exponents for one-dimensional cellular automata[END_REF], which he defines as a shift-invariant measure of left and right perturbation speeds. In fact he proves an inequality connecting the measure-theoretical entropy of a cellular automaton and its shift-invariant Lyapunov exponents. Tisseur altered the definitions slightly and considers average Lyapunov exponents in [START_REF] Tisseur | Cellular automata and lyapunov exponents[END_REF], where he establishes a similar connection. Finally the pointwise Lyapunov exponents are defined in [START_REF] Bressaud | On a zero speed sensitive cellular automaton[END_REF].

Connections between various dynamical properties of cellular automata and possible values of Lyapunov exponents have been studied for example in [START_REF] Finelli | Lyapunov exponents vs expansivity and sensitivity in cellular automata[END_REF], [START_REF] Manzini | On computing the entropy of cellular automata[END_REF] and [START_REF] Bressaud | On a zero speed sensitive cellular automaton[END_REF]. In [START_REF] Manzini | On computing the entropy of cellular automata[END_REF], a closed formula for calculating the value of the shift-invariant Lyapunov exponents, was presented in the setting of linear cellular automata. In [START_REF] Finelli | Lyapunov exponents vs expansivity and sensitivity in cellular automata[END_REF], the authors proved that for a given positively expansive cellular automaton, the shift-invariant Lyapunov exponents are positive for each configuration. The result was later improved in [START_REF] Bressaud | On a zero speed sensitive cellular automaton[END_REF], where it was shown that the same result holds for the pointwise Lyapunov exponents.

In [START_REF] Bressaud | On a zero speed sensitive cellular automaton[END_REF], Bressaud and Tisseur construct a sensitive cellular automaton, for which the value of the left and right average Lyapunov exponents is zero, with respect to a specific measure. In the same paper the Conjecture 3 states, that for a given sensitive cellular automaton it is necessary, that there exists such a configuration, whose either left or right pointwise Lyapunov exponent has a positive value. The same conjecture is reinstated by Kůrka as a Conjecture 11 in [START_REF] Kůrka | Topological Dynamics of Cellular Automata[END_REF]. Our aim is to prove this conjecture false, which we will do in Theorem 2. The result follows from the existence of complete, aperiodic Turing machines, established in [START_REF] Blondel | On the presence of periodic configurations in turing machines and in counter machines[END_REF] and [START_REF] Cassaigne | A small minimal aperiodic reversible turing machine[END_REF], and their movement bounds proved in [START_REF] Jeandel | Computability of the entropy of one-tape turing machines[END_REF] and [START_REF] Guillon | Distortion in one-head machines and cellular automata[END_REF].

Preliminaries

An alphabet Σ is a finite set of symbols. A word of length n over an alphabet Σ is any element w = (w 0 , w 1 , . . . w n-1 ) = w 0 w 1 • • • w n-1 from the set Σ [0,n) = Σ n and |w| = n is the length of a word w. The empty word is denoted as and it is the unique word of length zero. A set of all finite words i.e. n∈N Σ n is denoted as Σ * and a set of all finite non-empty words Σ * \{ } is denoted as

Σ + . A concatenation • : (Σ * ) 2 → Σ * is a mapping such that u • v = u 0 u 1 . . . u n v 0 v 1 . . . v m , where u = u 0 u 1 . . . u n and v = v 0 v 1 . . . v m .
We will adapt the shorthand notation uv for the concatenation of any two words. Elements from the sets Σ N , Σ Z- and Σ Z are called right-infinite, left-infinite and bi-infinite words, respectively. Furthermore we define a set

Σ Ω = Σ + ∪ Σ N ∪ Σ Z-∪ Σ Z . A concatenation of elements u ∈ Σ Ω and v ∈ Σ Ω is
defined when u is finite or left-infinite and v is finite or right-infinite. Let u ∈ Σ Ω and w ∈ Σ Ω , we will denote u w if there exists such j ∈ Z, that u i+j = w i for each i in the domain of u, and say that u is a subword of w. If Σ and Γ are two alphabets, we will denote the set

{uv | u ∈ Σ α , v ∈ Γ β } as Σ α Γ β .
where α ∈ {Z -, * , +}∪N and β ∈ {N, * , +}∪N. In this notation, if Σ = {a}, we will omit the brackets. Finally if w ∈ Σ * , we will use the notation w ∞ for the right-infinite word ww • • • .

A Turing machine is a 3-tuple (Q, Γ, δ), where Q is a finite set of states, Γ is a finite set of symbols and δ is a partial mapping δ :

Q × Γ → Q × Γ × ∆ called a
transition rule, where ∆ = {←, →}. A configuration is a 3-tuple (w, i, q), where w ∈ Γ Z , i ∈ Z and q ∈ Q. We will write (w, i, q) (w , j, r) if δ(q, w i ) = (r, w i , d), where w ∈ Γ Z , and j = i+1 if d = → and j = i-1 if d = ← and w k = w k for each k = i. Inductively we define n , where is applied n times. Furthermore we will write (w, i, q) + (w , j, r) if there exists such n ∈ Z + , that (w, i, q) n (w , j, r) holds. We will call a pair (q, a) ∈ Q × Γ an error pair if δ(q, a) is undefined. A Turing machine is called complete if it has no error pairs. We will call a configuration (w, i, q) periodic if (w, i, q) + (w, i, q) and weakly periodic if there exists such j ∈ Z, that (w, i, q) + (w , i+j, q), where w k+j = w k , for each k ∈ Z. We will call a Turing machine periodic if all its configurations are periodic, and aperiodic if none of its configurations are weakly periodic.

A topological dynamical system is a pair (X, f ), where X is a compact metric space and f is a continuous function f : X → X.

A shift dynamical system is a dynamical system (Σ Z , σ), where Σ is a finite set of symbols, Σ Z is the space called the full shift and σ, called the shift, is defined in a way that σ(x) i = x i+1 . The metric d of the space Σ Z is defined as d(x, y) = 2 -inf{|i|∈N|xi =yi} . It is not difficult to see that the space Σ Z is compact and that the function σ is continuous. An endomorphism is a continuous function

f : Σ Z → Σ Z , such that f • σ = σ • f . A one-dimensional cellular automaton is a 3-tuple A = (Σ, N, h), where Σ is a finite set of symbols called states, N is a neighbourhood (i 1 , i 2 , . . . , i n ) ∈ Z n and h : Σ n → Σ is a local rule. If N = [-r, r],
we call N a radius-r neighbourhood. In the context of cellular automata, we call the full shift Σ Z a configuration space and refer to its elements as configurations. The local rule together with the neighbourhood induces a global rule f : Σ Z → Σ Z , which is defined in such a way that f (c) i = h(c i+i1 , c i+i2 , • • • , c i+in ). We make no distinction between a cellular automaton and its global rule. By the Curtis-Hedlund-Lyndon theorem, the cellular automata, sometimes abbreviated as CA, are exactly the endomorphisms of the shift dynamical systems.

Definition 1. A dynamical system (X, f ) is sensitive if ∃ > 0 : ∀δ > 0 : ∀x ∈ X : ∃y ∈ B δ (x) : ∃n ∈ N : f n (y) ∈ B (f n (x)).
Definition 2. Let (Σ, N, h) be a one-dimensional cellular automaton, with a global rule f : Σ Z → Σ Z . For every c ∈ Σ Z , we define

W + m (c) = {c ∈ Σ Z | ∀i ≥ m, c i = c i } and W - m (c) = {c ∈ Σ Z | ∀i ≤ m, c i = c i }. Furthermore we define I + n (c) = min{m ∈ N | f i (W + -m (c)) ⊆ W + 0 (f i (c)), ∀i ≤ n} and I - n (c) = min{m ∈ N | f i (W - m (c)) ⊆ W - 0 (f i (c)), ∀i ≤ n}
Finally we define the pointwise Lyapunov exponents as

λ + (c) = lim inf n→∞ I + n (c) n and λ -(c) = lim inf n→∞ I - n (c) n .
Definition 3. Let X be a set. A relation is a subset R ⊆ X × X. We will use the standard notation aRb if (a, b) ∈ R. We will denote the complement of R as R c , i.e. R c = (X × X) \ R.

One way to simulate Turing machines in the CA setting, is to use the moving head model, or TMH for short, as introduced in [START_REF] Kůrka | On topological dynamics of turing machines[END_REF] and by adding arrows, which separate different simulation areas from each other. We will be using a slight variation.

In the following definition we will describe a function induced from a given Turing machine. The function can then be used to define different types of local rules for CA. The construction is using the conveyor belt model, which has been used previously at least in [START_REF] Guillon | Distortion in one-head machines and cellular automata[END_REF], albeit the notations might differ. Definition 4. Let M = (Q, Γ, δ) be a Turing machine. We will define an alphabet for the conveyor belt model. Elements from Q 2 will be called Turing machine heads and elements from T 2 will be called tape symbols. The set ∆ = {→, ←}, consisting of direction symbols, is used for locating the Turing machine heads. The conveyor belt model uses two layers of tape, where the second one is reversed and the two tapes are connected at the ends of the simulation areas. We will define a relation, which allows us to define these simulation areas as simulation words, each consisting of at most one Turing machine head and where directions are forced to point to the unique Turing machine head. Let

Γ 2 = Q 2 ∪ T 2 , where Q 2 = Γ 2 × Q ∪ Q × Γ 2 , T 2 = Γ 2 × ∆.
Define a relation R 2 in a following way: Let a ∈ Γ 2 and b ∈ Γ 2 , then

aR 2 b if      a ∈ Γ 2 × {→} ∧ b ∈ (Γ 2 × {→}) ∪ Q 2 ∨ a ∈ Q 2 ∧ b ∈ Γ 2 × {←} ∨ a ∈ Γ 2 × {←} ∧ b ∈ Γ 2 × {←}.

Now we can define a set of simulation words:

Σ M,2 = {w ∈ Γ Ω 2 | w j R 2 w j+1 , ∀j}.
In the following, in the case that the input word contains a Turing machine head, we will assume uvwxy to be such a decomposition of the input word, that w ∈ Q 2 and that the following conditionals hold:

If |uv| > 0, then |v| = 1. If |xy| > 0, then |x| = 1.
This decomposition is then unique. In the cases when the input word consists of only tape symbols the decomposition can be arbitrary.

From the transition rule of the Turing machine, we can construct a mapping m 2 : Σ M,2 → Σ M,2 in the following way:

m 2 (uvwxy) =                                                                        uv w xy if w = (q, a, a ) ∧ v = (b, b , →) ∧ δ(q, a) = (r, c, ←), where v = (r, b, b ) ∧ w = (c, a , ←), w xy if uv = ∧, w = (q, a, a ) ∧ δ(q, a) = (r, c, ←), where w = (c, a , r), uv w xy if w = (a, a , q) ∧ v = (b, b , →) ∧ δ(q, a ) = (r, c , →), where v = (b, b , r) ∧ w = (a, c , ←), w xy if uv = ∧, w = (a, a , q) ∧ δ(q, a ) = (r, c , →), where w = (r, a, c ), uvw x y if w = (q, a, a ) ∧ x = (b, b , ←) ∧ δ(q, a) = (r, c, →), where x = (r, b, b ) ∧ w = (c, a , →), uvw if xy = ∧, w = (q, a, a ) ∧ δ(q, a) = (r, c, →), where w = (c, a , r), uvw x y if w = (a, a , q) ∧ x = (b, b , →) ∧ δ(q, a ) = (r, c , ←), where x = (b, b , r) ∧ w = (a, c , →), uvw if xy = ∧, w = (a, a , q) ∧ δ(q, a ) = (r, c , ←),
where w = (r, a, c ), uvwxy otherwise.

Definition 5. Let M = (Q, Γ, δ) be a Turing machine. Using the notations introduced in Definition 4, we define a tracking function

T : Σ M,2 × N → Z, where T (w, j) = k if m j 2 (w) k ∈ Q 2 and T (w, j) = 0 if w i ∈ Q 2
for each i in the inputs domain. We also define the number of indices visited by a Turing machine head in j steps, given some initial word, by a function V (w, j) : Σ M,2 × N → N, defined as V (w, j) = |{T (w, j ) ∈ Z | j ≤ j}|. Finally we define a movement bound M : N → N as a function such that M (j) = max w∈Σ M,2 V (w, j).

The sublinearily of the movement bound for aperiodic Turing machines is already proved in [START_REF] Jeandel | Computability of the entropy of one-tape turing machines[END_REF]. We will however make use of the following tighter bound. Theorem 1. [START_REF] Guillon | Distortion in one-head machines and cellular automata[END_REF] Let M = (Q, Γ, δ) be Turing machine, and M its movement bound. If M is aperiodic, then M = O( n log n ).

Lyapunov exponents for sensitive cellular automata

In this section we study the question whether or not the sum of the pointwise Lyapunov exponents of a given sensitive CA can take the value 0 or not. We use the notations of the previous sections without explicitly referring to them.

In Theorem 2 we show how to construct a sensitive cellular automaton for a given aperiodic Turing machine, such that its left and right pointwise Lyapunov exponents are bounded from above by a sublinear function derived from the movement bound of the Turing machine.

Theorem 2. There exists such a sensitive one-dimensional cellular automaton

(Σ, N, f ), that λ + (c) = λ -(c) = 0 for every configuration c ∈ Σ Z .
Proof. Let M = (Q, Γ, δ) be an aperiodic, complete Turing machine. We will alter M slightly to produce a new Turing machine M = (Q , Γ, δ ), where Q = Q × {0, 1} and δ ((q, x), a) = ((r, x), b, d), where δ(q, a) = (r, b, d). Essentially this modification gives us two copies of the original Turing machine and the computations of the new machine can be projected to the computations of the original machine. Using the new Turing machine, we can set our symbol set for the CA, which we are constructing, as Σ = Γ 2 ∪ {>}. We will refer the symbol > as the eraser. We will also use the notations

Q i 2 = Γ 2 × Q × {i} ∪ Q × {i} × Γ 2 , where i ∈ {0, 1} and Q 2 = Q 0 2 ∪ Q 1 2 .
For each configuration c ∈ Σ Z , we define a set of locations for the Turing machine heads as

H c = {i ∈ Z | c i ∈ Q 2 }
and a set of locations for the erasers as

E c = {i ∈ Z | c i = >}.
Next we define the simulation bounds as functions l c : H c → Z ∪ {-∞} and r c : H c → Z ∪ {∞} in the following way:

l c (i) = sup{j ∈ Z | j ≤ i and c j-1 R c 2 c j } and r c (i) = inf{j ∈ Z | i ≤ j and c j R c 2 c j+1 }.
From these bounds we can define the set of cells that are not part of any simulation area as

U c = Z \ ( i∈Hc (l c (i) -1, r c (i) + 1)).
Using the simulation bounds, we can define a function, which simulates all Turing machines in their respective simulation areas as m : Σ Z → Σ Z , where

m(c) (lc(i)-1,rc(i)+1) = m 2 (c (lc(i)-1,rc(i)+1) ) ∀i ∈ H c and m(c) i = c i ∀i ∈ U c .
Clearly m is a cellular automaton since we can extract a radius-1 local rule from its definition.

Next we fix a state a 0 ∈ Γ 2 × {←} and define a second cellular automaton e : Σ Z → Σ Z , whose local rule is described in Table 1. The dynamics of e can be described in the following way: If a cell i, containing the eraser state >, sees a Turing machine head from the state set Q 1 2 at a cell i -1, then the eraser state gets pushed to the next cell at i + 1. This happens regardless of the previous content of the cell i + 1. As a consequence, the simulation area in the left side increases by one cell. If there is a simulation area on the right side, then either its size is reduced by one or it is removed entirely. The simulation area gets removed if there is a Turing machine head at a cell i + 1 or i + 2. In such situation, the Turing machine head is replaced with a tape symbol a 0 . Otherwise if a cell at state > does not see a Turing machine head from the state set Q 1 2 to their immediate left, then it stays at its current state. When a Turing machine head from the state set Q 1 2 moves an eraser state, it changes to an equivalent symbol from the state set Q 0 2 . If a Turing machine head from the state set Q 0 2 visits the left bound of the simulation area, which it belongs to, then if it does not get erased, it changes to an equivalent symbol from the state set Q 1 2 . The idea is that the simulation areas can increase in their size arbitrarily far to the right, but only one cell at a time and in between the increments, the Turing machine head must visit the left bound of the simulation area.

q, q ∈ Q2, qi = (q, i), x ∈ Σ u ∈ Σ \ Q 1 2 , u ∈ Σ \ Q2, v, v ∈ Σ \ {>}, v ∈ Γ 2 × {←} ∪ Q2 q1 > a0 x > u u q1 > a0 x > q a0 v v q1 q0 > a0 u > > u u u > > q a0
v v q0 q1 Table 1. The definition of the local rule of the CA e : Σ Z → Σ Z . Here a0 is a fixed state from the set Γ 2 × {←}. The rule uses a neighbourhood (-2, -1, 0, 1). Input is written on the first row and output on the second row. If the output is written, but the input is partially missing, it means that the missing cells do not affect the output. For the inputs that do not appear in the table, the local rule behaves as the identity mapping.

We are ready to define our cellular automaton of interest as f : Σ Z → Σ Z , where f = m • e. The behaviour of the CA is depicted in Figure 1.

In Lemmas 1 and 2 we prove, that the CA we constructed has the desired properties and therefore the claim follows.

Lemma 1. The cellular automaton (Σ, N, f ) constructed in Theorem 2 has the property that λ + (c) = λ -(c) = 0 for every configuration c ∈ Σ Z .

Proof. We will begin the proof by introducing tracking functions for the eraser states, the simulation bounds and the Turing machine heads. The point of them is, that given an initial configuration and a cell containing a Turing machine head, an eraser or a simulation bound, we can tell to which cell said symbol has travelled to in time. In the output of each function, we will use the symbol -to denote, that the symbol no longer exists, i.e. it has been destroyed by a symbol >. We will then show that a difference can only propagate inside each simulation area or by the movement of the simulation areas. In either case the movement is bounded from above by a sublinear function derived from the movement bound of the Turing machine. 2 and Q 0 2 , respectively, and the red lines represent the erasers. We also note, that in the figure, time increases from top to bottom. One can witness several types of behaviour in the two simulation areas. When the Turing machine head in the state from the set Q 1 2 , of the left simulation area, visits the right boundary, then in the next time-step, the eraser moves one cell to the right, which also moves the left boundary of the second simulation area. In the same time step, the Turing machine head switches to some state in the set Q 0 2 . Such Turing machine heads do not move the eraser states as witnessed when the cyan coloured line visits the right boundary. For the Turing machine to be allowed to move the eraser state again, it needs to switch back to an element from Q 1 2 , which happens if and only if it visits the left boundary. The right simulation area does not have an eraser on the right side and hence its simulation area can never increase in size. In the middle of the image we can see that the Turing machine head on the right simulation area visits a cell within the distance of two, to a cell containing an eraser state and hence gets removed. This happens eventually in all simulation areas, in which at some time-step an eraser state can be seen in a cell within distance two to its left boundary.

First we will define the tracking function for the eraser states as e

c : E c ×N → Z ∪ {-}, where e c (i, j) =          i if j = 0, e c (i, j -1) + 1 if f j-1 (c) [ec(i,j-1)-1,ec(i,j-1)] ∈ Q 1 2 >, - if f j-1 (c) [ec(i,j-1)-2,ec(i,j-1)-1] ∈ Q 1 2 >, e c (i, j -1)
otherwise.

Next we will define the tracking functions for the simulation bounds and for the Turing machine head, inductively with respect to j ∈ N, as l c :

H c × N → Z ∪ {-}, r c : H c × N → Z ∪ {-} and h c : H c × N → Z ∪ {-}, where
and l c (i + 1) = i + 1 or l c (i + 1) = i and r c (i -1) = i -1. We only prove the case of the former as the latter is analogous. Immediately it follows that r c (i -1, 1) = i = r f (c) (i ) as i + 1 = >. The analysis that results in to showing that l c (i -1, 1) = l f (c) (i ) is similar to the case where there was only one Turing machine head in the interval [i -1, i + 1]. Hence we have that h c (i, 1) = i , where i = i -1. The case when i -1 ∈ H c , i ∈ H c , but i + 1 ∈ H c and the case when i -1 ∈ H c , i ∈ H c , but i + 1 ∈ H c are proved similarly. Now again by a straightforward induction it can be proved that for each n ∈ N and i ∈ H

f n (c) there exists unique i ∈ H c , such that l c (i, n) = l f n (c) (i ), r c (i, n) = r f n (c) (i ) and h c (i, n) = i .
Denote by S c (i, j) = (l c (i, j) -1, r c (i, j) + 1), w i,j = c Sc(i,j) and s i,j = |w i,j |. Suppose that i ∈ H c and j ∈ N are such that r c (i, j + 1) = r c (i, j) + 1. By definition of the tracking function, we have that f j (c) [rc(i,j),rc(i,j)+1] ∈ Q 1 2 > and by the definition of the CA e, we have that f j+1 (c) hc(i,j+1) ∈ Q 0 2 . We assume that there exists such k ∈ N, that r c (i, j + 1) < r c (i, j + 1 + k). Then necessarily, from the aperiodicity of the Turing machine, we get that there exists a minimal such k ≥ 1, that h c (i, j + k ) = l c (i, j + k ). From the definition of e, this means that f j+k (c) hc(i,j+k

) ∈ Q 0 2 for each 1 ≤ k ≤ k and f j+k +1 (c) hc(i,j+k +1) ∈ Q 1 2 .
Notice that it also implicitly holds that > f t (c) [lc(i,t)-2,lc(i,t)-1] for each t ≤ j + k as otherwise we would have h c (i, j + k ) = -. We denote by

M -1 (n) = min{t ∈ N | M (t) ≥ n} for each n ∈ N.
By the movement bound of the Turing machine we have that M -1 (s i,j ) ≤ k ≤ p si,j+1 , where p = |Σ|. Furthermore from the definition of k, we have that M -1 (s i,j+1 ) ≤ k -k ≤ p si,j+1 and hence 2M -1 (s i,j ) ≤ k ≤ 2p si,j+1 . It is easy to see that during these k steps a valid Turing machine computation is performed with the input word w i,j+k . Finally let z j = s i,j -s i,0 and assume that r c (i, j) = h c (i, j). From the above considerations we get that M -1 (s i,0 + z j -1) ≤ j ≤ 2 1≤k≤zj p (si,0+k) ≤ p (si,0+2+zj ) . From the lower bound it follows that |{α i,j ∈ N | j ≤ n}| ≤ M (n), for each n ∈ N and i ∈ H c ∪ E c , where α is any of the tracking functions for either the simulation bounds, the Turing machine heads or the erasers.

We are ready to prove the claim of the theorem. Let c ∈ Σ Z . For an infinite set of integers n ∈ N, we want to find such values m n ∈ N, that if c ∈ W + -mn (c) and there exists j n , such that f jn (c ) ∈ W + -mn (f jn (c)), then j n > n. For each n, we can first assume that m n ≥ M (n). We will analyse three cases: 1) None of the cells in the interval [-m n -3, 0] are in the eraser state in either of the configurations, during the first n iterations. 2) There is an upper bound k < n, after which none of the cells in the interval [-m n -3, 0] are in the eraser state in either of the configurations. 3) One of the configurations has a cell in the interval [-m n -3, 0], which is at an eraser state at the n th iteration.

First let us assume that > f j (c ) [-mn-3,0] holds for each j ≤ n and c ∈ {c, c }. Then by the definition of the global rule e, we have that e(f j (c )

) i = f j (c ) i for each j ≤ n, i ∈ [-m n , 1] and c ∈ {c, c }. Hence if j n ≤ n, is such that f jn (c ) ∈ W +
-mn (f jn (c)), then the difference must be caused by applying the global rule m. By the definition of m, this means that f jn (c ) 0 ∈ Q 2 for either c = c or c = c. From the movement bound and the assumption that m n ≥ M (n), it then follows that there exists such i ∈ H c , that h c (i, j) ∈ [-m n , m n ] for each j ≤ n. Since c ∈ W + -mn (c), then i ∈ H c ∩ H c and hence locally within the cells in [-m n , m n ], the same Turing machine computation is simulated for the n steps and thus j n > n.

Assume then that there exists such j ≤ n and c ∈ {c, c }, that > f j (c ) [-mn-3,0] . We further assume, that there exists minimal such k < n, that > f j (c ) [-mn-3,0] for each c ∈ {c, c } and j ∈ N, such that k < j ≤ n. Then due to definition of e, necessarily for either

c = c or c = c it holds that f k (c )[-1, 0] ∈ Q 1 2 > and > f k (c ) [-mn-3,1]
. We can assume that this holds for c = c. Then, from the movement bound it follows, that there ex-

ists such i ∈ H c , that h c (i, j) ∈ [-m n -1, -1] for each j ≤ k. We clearly have that i ∈ H c . If h c (i, j) = h c (i, j)
for each j ≤ n, then a difference cannot propagate to the origin during the n steps as again the same computation would happen in both configurations. On the other hand the only way that h c (i, j) = h c (i, j) for some j ≤ n is if h c (i, j) = -. This is however impossible as it would require that > c [-mn-3,i-1] , which would mean that > f k+1 (c ) [-mn-3,0] , which is against the assumption. Therefore again if j n ∈ N, is such that f jn (c ) ∈ W + -mn (f jn (c)), then from the above consideration, we have that j n ≥ n.

Finally, for the last case we assume that > f n (c) [-mn-3,0] . Then from the movement bound M , it follows that there exists such i e ∈ E c , that e c (i e , j) ∈ [-2M (n) -3, 0] for each j ≤ n. We can then assume that m n ≥ 2M (n) + 3 and hence i e ∈ E c . If e c (i e , j) = e c (i e , j) for each j ≤ n, then directly from the definition of the CA f , it follows that f j (c) k = f j (c ) k , for each j ≤ n and k ≥ e c (i e , j) and hence in such case j n > n. Therefore we consider the case where there exists such j ≤ n, that e c (i e , j) = e c (i e , j). Suppose that for each pair i ∈ H c and i ∈ E c , such that i < i = r c (i) + 1 < i e , it holds that either there exists minimal such j 0 ≤ n, that e c (i , j 0 ) -2 ≥ h c (i, j 0 ) and e c (i e , j) = e c (i e , j 0 ) for each j 0 < j ≤ n or e c (i , j) -2 < h c (i, j) for each j ≤ n, where c ∈ Σ Z is such a configuration that c j = c j for each j > l c (i) -1 and c j = a 0 for each j < l c (i). If we now assume that m n ≥ 3M (n) + 4, then any i ∈ E c , such that i < -m n cannot effect the value of e c (i e , j) for any j ≤ n and neither can any i ∈ E c ∩ [-m n , i e -1] by our assumption regarding configurations c . Hence if e c (i e , j) = e c (i e , j) and j ≤ n, the difference is caused by a Turing machine head. But during the first n steps said Turing machine head can only visit cells [-M (n

) + i e -1, i e + M (n)] ⊆ [-3M (n) -4, M (n)]
. So again both configurations c and c are simulating the exact same computation during the first n steps and thus j n > n.

Let us assume that there exists such k ∈ N and indices i j ∈ H c , where 1 ≤ j ≤ k + 1, that i j < i j+1 for each j ≤ k. Let us denote as c j ∈ Σ Z such a configuration that c j j = c j for each j > l c (i j ) -1 and c j j = a 0 for each j < l c (i j ). We will also assume that for each j ≤ k + 1, we have that e ij = r c (i j ) + 1 ∈ E c and we will further assume that there exists such an increasing finite sequence of times t j , that h c j (i j+1 , t j ) ∈ [e c j (e ij , t j ) + 1, e c j (e ij , t j ) + 2] and e c j+1 (e ij+1 , t j+1 ) > e c j+1 (e ij+1 , t j ) for each j ≤ k. That is the leftmost simulation area in c j , destroys the simulation area of the Turing machine head i j+1 at time-step t j . In c however, this might not happen as there could be another simulation, which destroys the leftmost simulation area of c j before time t j . This could allow an existence for an alternating chain of such simulation areas, where simulation areas of i j are destroyed for each even j in configuration c and and odd j in configuration c . We want to find an upper bound a(n) for the maximum distance e c k (e i k , t i k ) -e c 1 (e i1 , t i1 ), assuming that t k ≤ n, as then we know how fast a difference can potentially propagate via such a chain of simulation areas. In such a case, suppose that we would have that c

[ie-a(n)-2M (n)-4,∞) = c [ie-a(n)-2M (n)-4,∞)
, where i e = e i k+1 . Assume that there exist such j 0 ∈ N, that e c (i e , j 0 ) = e c (i e , j 0 ), then if the difference is due to a chain that we have described above, it must be that in some of the simulation areas of such chain, a different computation is performed in the configurations c and c . If i e -e c k (e i k , t i k ) > M (n) + 2, then we would have that e c (i e , j) = e c (i e , j) for each j ≤ n. If i e -e c k (e i k , t i k ) ≤ M (n) + 2, then e c 1 (e i1 , t i1 ) ≥ i e -2-M (n)-a(n), but then in all of the simulation areas the same computation is performed during n iterations and hence e c (i e , j) = e c (i e , j) for each j ≤ n.

Let us assume that for each j ≤ k, we have that p n ≤ t j ≤ p n+1 , where p = |Σ| and n ≥ 5. Let b j = e c j+1 (e ij+1 , t j+1 )-e c j (e ij , t j ) and b j,t = e c (e ij+1 , t)e c (e ij , t). Since h c j (i j+1 , t j ) ∈ [e c j (e ij , t j )+1, e c j (e ij , t j )+2] and e c j+1 (e ij+1 , t j+1 ) > e c j+1 (e ij+1 , t j ) for each j ≤ k, there must exist such j ≤ t j+1 -t j , that h c j+1 (i j+1 , t j + j ) = e c j+1 (e ij+1 , t j+1 ) -1 and therefore we have that b j -2 ≤ M (t j+1 -t j ) for each j ≤ k.

By Theorem 1 there exists such h : N → N, that M (n) ≤ h(n) = C n log(n) for each n ∈ N and where C > 0. It is easy to see that there exists such a positive real number C , that h is concave in the domain [C , ∞). We will split the times t j into two sets A and B, such that j ∈ A if t j+1 -t j ≥ C and t j ∈ B otherwise. We have that

k j=1 b j ≤ 2k + k j=1 M (t j+1 -t j ) ≤ 2k + k j=1 h(t j+1 -t j ) ≤ 2k + |B|h(C ) + tj ∈A h(t j+1 -t j ) ≤ k(2 + h(C )) + |A|h( t j ∈A tj+1-tj |A| ) ≤ k(2 + h(C )) + |A|h( p n+1 |A| ) = k(2 + h(C )) + C p n+1 log( p n+1 |A| ) ≤ k(2 + h(C )) + C p n+1 log( p n+1 k )
,

where the the fourth inequality follows from Jensen's inequality for concave functions. Hence the upper bound is maximized when k is maximized. Thus we want to find an upper bound for the value k. First for each t < t 1 , we have that where b 0 = e c (e i1 , t 1 ) -e c (e i1 , t). From the movement bounds we proved earlier and since for each j we assumed that t j ≥ p n , we have that b j,p n-1 > n -3 for each j ≤ k. This follows from the fact that the Turing machine head i j is not erased before p n steps in the configuration c j-1 and p n-1 is enough time to have have had any Turing machine head visit n -3 cells, even if the simulation area had started from a size 1. Therefore it follows that ≤ 2C p n+1 log(n-3) .

Recall that t j were assumed to be inside the interval [p n , p n+1 ], where n ≥ 5. For the times less than p 5 , we get some constant upper bound C and hence taking union of intervals [p i , p i+1 ], where 1 ≤ i ≤ n, we get our upper bound

a(n) = C + n i=5 C p i+1 log(i) ≤ C + C p n+2 log(n + 1)
.

Hence if we choose m n = 5M (p n+1 ) + C + C p n+2 log(n+1) , we have that e c (i e , j) = e c (i e , j) for each j ≤ p n+1 .

Combining all the three cases, we have shown that Let us then suppose that I = {i ∈ N | f i (c) 0 = >} is infinite. Then for each k < 0, we choose c ∈ W + k (c), such that c [k-1,k] ∈ Q 1 2 > and c i = a 0 for each i < k -1. From the aperiodicity it again follows that there exists such n ∈ N that > f n (c ) (-∞,0] for each n > n.

Fig. 1 .

 1 Fig.1. In this figure, we have depicted the behaviour of the CA constructed in Theorem 2. The black lines represent the left and right simulation bounds, the blue and cyan lines represent Turing machine heads from the sets Q 1 2 and Q 0 2 , respectively, and the red lines represent the erasers. We also note, that in the figure, time increases from top to bottom. One can witness several types of behaviour in the two simulation areas. When the Turing machine head in the state from the set Q 1 2 , of the left simulation area, visits the right boundary, then in the next time-step, the eraser moves one cell to the right, which also moves the left boundary of the second simulation area. In the same time step, the Turing machine head switches to some state in the set Q 0 2 . Such Turing machine heads do not move the eraser states as witnessed when the cyan coloured line visits the right boundary. For the Turing machine to be allowed to move the eraser state again, it needs to switch back to an element from Q 1 2 , which happens if and only if it visits the left boundary. The right simulation area does not have an eraser on the right side and hence its simulation area can never increase in size. In the middle of the image we can see that the Turing machine head on the right simulation area visits a cell within the distance of two, to a cell containing an eraser state and hence gets removed. This happens eventually in all simulation areas, in which at some time-step an eraser state can be seen in a cell within distance two to its left boundary.

  t = e c (e i k+1 , t) -e c (e i1 , t) ≤ e c (e i k+1 , t k+1 ) -e c (e i1 , t) = e c (e i k+1 , t k+1 ) -e c (e i1 , t 1 ) + e c (e i1 , t 1 ) -e c (e i1 , t) ≤ e c (e i1 , t 1 ) -e c (e i1 , t) +

bM - 1

 1 j ≥ k(n -3). The function M -1 that we introduced earlier is bounded from below by identity mapping. Therefore we have thatp n+1 ≥ k j=0 (b j ) = k(n -3).Therefore we have the upper bound k ≤ p n+1 n-3 . Denoting C = max{C, h(C )+2} and by combining our inequalities we have that k j=0 b j ≤ k(2 + h(C )) + C p n+1 log( p n+1 k ) ≤ C p n+1 n-3 + C p n+1 log(n-3)

Lemma 2 .

 2 0 as n goes to infinity. As this holds for each c ∈ Σ Z and n ∈ N. Hence we have that λ + (c) = 0 for each c ∈ Σ Z .The fact that λ -(c) = 0 holds for each c ∈ Σ Z is much easier to see. First of all we have seen that the eraser states travel only to the right direction. Hence if there exists such i ∈ N thatc i = >, for some i ≥ 0, it means that if c ∈ W - i (c), then f n (c ) ∈ W - i (f n (c)),for each n ∈ N. Hence any difference coming from right must propagate within a single simulation area. But then it follows from the movement bound that I - n (c) ≤ M (n) and hence λ -(c) = 0 for each c ∈ Σ Z . The cellular automaton (Σ, N, f ) constructed in Theorem 2 is sensitive.Proof. Let c ∈ Σ Z and suppose thatI = {i ∈ N | f i (c) 0 = >} is finite. Then there exists such m ∈ N, that f n (c) 0 = > for each n ≥ m. Let k ≤ -M (m) -1 and c ∈ W + k (c), such that c [k-1,k] ∈ Q 1 2 > and c i = a 0 for each i < k -1.We saw in the proof of Lemma 1, that e c (k, p n+2 ) -e c (k, 0) ≥ n for each n ∈ N and hence there exists such n ≥ m, that e c (k, n) = 0.

Acknowledgements

The author acknowledges the emmy.network foundation under the aegis of the Fondation de Luxembourg for its financial support.

if j = 0, l c (i, j -1) + 1 if f j-1 (c) [lc(i,j-1)-2,lc(i,j-1)+1] ∈ Q 1 2 > T 2 T 2 , l c (i, j -1) if f j-1 (c) [lc(i,j-1)-2,lc(i,j-1)-1] ∈ Q 1 2 >, and

We will show that for each i ∈ E f (c) there exists unique i ∈ E c , such that e c (i, 1) = i . By the definition of f it follows immediately, that c i = > or c i -1 = > and furthermore if only the latter holds, then necessarily

2 >, then i = i + 1 and otherwise i = i . Then furthermore by a straightforward induction, it holds that for each n ∈ N and i ∈ E f n (c) there exists unique i ∈ E c , such that e c (i, n) = i .

Analogously we will show that for each i ∈ H f (c) there exists unique i ∈ H c , such that h c (i, 1) = i , l c (i, 1) = l f (c) (i ) and r c (i, 1) = r f (c) (i ). From the definition of the CA m, it follows immediately that there exists such i ∈

But this is consistent with the definition of l c (i, 1) and hence l c (i, 1) = l f (c) (i ). Similarly we can show that r c (i, 1) = r f (c) (i ). It is then apparent, that h c (i, 1) = i as cell i is the only cell containing a Turing machine head in the interval (l c (i, 1) -1, r c (i, 1) + 1). If there exist multiple Turing machine heads in the interval [i -1, i + 1], then let us first assume that i -1 ∈ H c and i + 1 ∈ H c . Then if i ∈ H c it follows that l c (i ) = i = r c (i ) and since i -1 = > and i + 1 = >, we have that l c (i , 1) = i = r c (i , 1) and hence the claim follows. If i ∈ H c , then i ∈ T 2 , hence either we have that r c (i -1) = i


Gabriel Peyré

DYNAMIC TEXTURE SYNTHESIS WITH GROUPLETS

This paper proposes a new method to synthesize dynamic geometric textures from a given exemplar image. It solves the Navier Stokes equations of fluid dynamics to drive the evolution of a geometric layer and a grouplet coefficients layer. The geometric layer describes the turbulent structures of the texture while the coefficient layer lays a set of strokes along the flows. These strokes corresponds to grouplet atoms, that have a wide range of width and lengths. These grouplets atoms generate a tight frame that is adapted to represent the texture to analyze and synthesize.

Résumé. Ce papier décrit une nouvelle méthode pour synthétiser des textures géométriques dynamiques à partir d'une image donnée en exemple. Les évolutions d'un flot géométrique et des coefficients de grouplets suivent les équation de Navier Stokes. Le flot géométrique dicte la structure turbulente de la texture, tandis que les coefficients de grouplets génèrent des formes le long de ce flot. Les atomes de grouplets ont un support allongé suivant une large gamme de longueurs et de largeurs. Ces atomes forment une trame ajustée qui est adaptée à la représentation compacte des textures à analyser et à synthétiser.

Natural images often contain regions composed of locally oriented structures, such as those depicted in Figure 1. These anisotropic textures are said to be locally parallel since they are composed of approximately parallel oscillations that propagate over the image plane. This paper is focussed on the animation of a static locally parallel texture according to fluid dynamics motions. This creates a turbulent texture whose patterns are similar to those of the exemplar.

Introduction

Anisotropic texture processing. The analysis of the local geometry of textures is studied extensively in computer vision through the computation of local differential estimators, see for instance [START_REF] Kass | Analyzing oriented patterns[END_REF][START_REF] Kothe | Edge and junction detection with an improved structure tensor[END_REF]. Section 2 uses such a local orientation descriptor to build the geometric layer of the grouplet representation.

Geometric decompositions. Texture models based on non-adaptive wavelet decompositions fail to represent in a compact manner geometric singularities. Geometrical decompositions improve over the wavelets for the approximation of edges and textures in images. Local oscillating atoms such as Gabor [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF], steerable wavelets [START_REF] Simoncelli | Shiftable multi-scale transforms[END_REF] or wave-atoms [START_REF] Demanet | Wave atoms and sparsity of oscillatory patterns[END_REF] better capture the directionality of textures. These fixed transforms are however not adaptive and are sub-optimal for texture processing. Adaptive representations such as bandlets [START_REF] Spencer | Bandelet Image Approximation and Compression[END_REF][START_REF] Mallat | Orthogonal bandlets for geometric image approximation[END_REF] and grouplets [START_REF] Mallat | Geometrical grouplets[END_REF] are tuned for a specific image to process. This adaptivity is achieved by computing an optimized geometry that parameterizes the representation. Section 2.2 reviews the grouplet transform that computes the coefficients layer of the grouplet representation.

Static texture synthesis. Texture synthesis is performed by sampling a texture model that constrains the geometric patterns of the images. Fourier and fractals modeling of texture [START_REF] Lewis | Texture synthesis for digital painting[END_REF][START_REF] Perlin | An image synthesizer[END_REF] creates cloud-like texture that are suitable to model some natural phenomena. Multiscale models constrain the distribution of wavelet coefficients of textures that characterizes point-wise singularities in textures. Retaining only wavelets marginal coefficients generates textures without geometric patterns [START_REF] Cook | Wavelet noise[END_REF][START_REF] Heeger | Pyramid-Based texture analysis/synthesis[END_REF]. Higher order statistical modeling of wavelets coefficients [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF][START_REF] Zhu | Filters, random fields and maximum entropy (FRAME): Towards a unified theory for texture modeling[END_REF] improves the visual quality of synthesis by capturing geometric singularities. The representation of local features in these methods is however implicit and thus hard to analyze or control.

Computer graphics methods generate new images of a given texture through careful and consistent copying of pixels and patches from an example image [START_REF] Efros | Texture synthesis by non-parametric sampling[END_REF][START_REF] Efros | Image quilting for texture synthesis and transfer[END_REF][START_REF] Wei | Fast texture synthesis using tree-structured vector quantization[END_REF]. Despite the high visual quality of the results, these approaches do not provide a parametric model for compact charactization of texture classes. Dynamic Textures. Fluid texture synthesis generates oriented patterns that are warped along a coherent flow. Such a turbulent flow is animated by discretizing the Navier-Stokes equations of fluid motion [START_REF] Foster | Realistic animation of liquids[END_REF][START_REF] Stam | Stable fluids[END_REF] or using texture synthesis of vector fields [START_REF] Bridson | Curl-noise for procedural fluid flow[END_REF]. Dynamic textures are rendered by advecting particules [START_REF] Neyret | Advected textures[END_REF], and is used to simulate physically plausible phenomena [START_REF] Lamorlette | Structural modeling of natural flames[END_REF][START_REF] Stam | Depicting Fire and Other Gaseous Phenomena Using Diffusion Processes[END_REF]. Wavelet domain texture synthesis enhances the visual quality of turbulence simulation by adding fine scale details [START_REF] Kim | Wavelet turbulence for fluid simulation[END_REF]. Dynamic textures can be generated by statistical modeling that captures the space and time homogeneity of natural dynamics [START_REF] Doretto | Dynamic textures[END_REF][START_REF] Kwatra | Graphcut textures: image and video synthesis using graph cuts[END_REF]. Patch recopy methods in computer graphics can be animated according to a vector field [START_REF] Kwatra | Texture optimization for example-based synthesis[END_REF][START_REF] Lefebvre | Appearance-space texture synthesis[END_REF] and can be mapped on a fluid in motion [START_REF] Bargteil | A texture synthesis method for liquid animations[END_REF][START_REF] Kwatra | Texturing fluids[END_REF]. Texture mixing animates an image by interpolating between two texture models [START_REF] Bar-Joseph | Texture mixing and texture movie synthesis using statistical learning[END_REF]. A grouplet model is introduced by Peyré [START_REF] Peyré | Texture synthesis with grouplets[END_REF] to perform geometric texture inpainting and synthesis. Section 4 uses the grouplet model to synthesize dynamic anisotropic textures.

Grouplet Texture Analysis

The grouplet framework is a two layers model composed of a geometric layer Γ and a coefficients layer that stores the decomposition of the texture in a waveletgrouplet tight frame B w (Γ). The wavelet-grouplet tight frame B w (Γ) = {b Γ j,ℓ,m } j,ℓ,m is composed of elongated atoms. Each b Γ j,ℓ,m is an anisotropic stroke of width ∼ 2 j and length ∼ 2 ℓ , located around a pixel m, and that follows the flow Γ.

The estimation of the model from an exemplar f 0 is performed by computing an optimized flow Γ adapted to f 0 as explained in Section 2.1, and the collection of coefficients { f 0 , b Γ j,ℓ,m } j,ℓ,m of f 0 in the wavelet-grouplet family B w (Γ), as explained in Section 2.3. The application of the model to perform dynamic texture synthesis corresponds to the generation of new geometric and coefficients layer, see Section 4.

Computation of the Geometric Layer

The geometric layer of our grouplet model is a vector field Γ the describes the local direction of a given input texture f 0 ∈ R N of N pixels. It is then discretized to obtained a set of discrete association field A ℓ that integrate the flow Γ on a distance of 2 ℓ . These discrete association fields drive the grouplet transform.

Local orientation and direction computation. The local direction of the edges in an image f 0 is captured by the vector v(x) orthogonal to the gradient ∇ x f 0 . This vector is computed numerically on the discrete grid {0, . . . , n -1} 2 using finite differences. While v is suitable to estimate the direction of step edges in images, it cannot be used directly to estimate the orientation of locally parallel textures such as those depicted on Figure 1. Indeed, the gradient vector vanishes at the top of ridges or at the bottom of valleys. As described in [START_REF] Kass | Analyzing oriented patterns[END_REF][START_REF] Kothe | Edge and junction detection with an improved structure tensor[END_REF], this problem is alleviated by averaging locally the tensor field vv T ,

T f0 (x) = (G σ ⋆ (vv T ))(x) = G σ ⋆ v 2 1 v 1 v 2 v 1 v 2 v 2 2 (x),
where G σ is a Gaussian kernel of variance σ 2 and where the convolution is applied on each component of the tensor. This scale σ should match the width of the texture oscillations, and is set to 5 pixels in the numerical experiments.

Each symmetric tensor T f0 (x) is decomposed as a sum of two rank-1 tensors

T f0 (x) = λ(x)(Γ 0 (x)Γ 0 (x) T ) + λ ⊥ (x)(Γ ⊥ 0 (x)Γ ⊥ 0 (x) T ). ( 1
)
where the eigenvalues are λ(x) λ ⊥ (x) 0 and (Γ 0 (x), Γ ⊥ 0 (x)) are the corresponding orthonormal eigenvectors. The dominant eigenvector Γ 0 (x) indicates the local direction of the texture.

This orientation field Γ 0 is not directional since Γ 0 (x) or -Γ 0 (x) can be used indifferently in the decomposition [START_REF] Bar-Joseph | Texture mixing and texture movie synthesis using statistical learning[END_REF]. A directional vector field Γ(x) = ε(x)Γ 0 (x) with ε(x) ∈ {-1, +1} is optimized to align the singularity of Γ with integral lines of the flow Γ 0 , see [START_REF] Peyré | Texture synthesis with grouplets[END_REF]. This geometrical flow Γ is the geometric layer that drives the computation of the grouplet coefficients layer.

Discrete association field. A set of non-local discrete association fields

{A ℓ } L-1 ℓ=0
is computed from the geometrical flow Γ. Each field A ℓ is then used to compute the grouplet decomposition at a scale 2 ℓ . Each association field A ℓ groups together two pixels x → y = A ℓ (x) that are separated by a distance of ||x -y|| ≈ 2 ℓ .

Our construction of the association fields A ℓ differs from the original one of [START_REF] Mallat | Geometrical grouplets[END_REF]. The original construction is based on a fixed directional ordering of the grid pixels, that forbids an arbitrary association field. Such a directional ordering is not well suited to process turbulent textures that might exhibit circular patterns or vortices.

The directional vector field Γ is the seed for a linear differential flow

∀ x, ∀ t > 0, dϕ x dt (t) = Γ(ϕ x (t)) and ϕ x (0) = x.
Each integral curve {ϕ x (t)} t>0 follows the vector field Γ. The discrete association field A ℓ at various scales ℓ is computed by sampling this integral curve at dyadic locations ∀ x, A ℓ (x) = [ϕ x (2 ℓ )] n×n ∈ {0, . . . , n -1} 2 , where [ • ] n×n is the rounding operator that projects the points on the sampling grid {0, . . . , n -1} 2 . Although the association fields are quantized in order to link points on the grid {0, . . . , n -1} 2 , they follow closely the structure of the input texture f 0 .

Grouplet Transform

Once the geometric layer Γ is known, the coefficients layer is computed with a fast redundant grouplet transform. A redundant grouplet transform is introduced by Mallat in [START_REF] Mallat | Geometrical grouplets[END_REF]. It corresponds to the decomposition of an image on a redundant family of vectors B(Γ). This family is a tight frame of R N chosen adaptively to process in an optimized way some given input image of N = n 2 pixels. This section describes a new construction of grouplets frame parameterized by a local geometric flow Γ. The grouplet atoms follow closely the flow lines of Γ. These grouplets are thus efficient to represent textures whose geometry is locally parallel to the orientation field Γ.

A grouplet tight frame

B(Γ) = b Γ ℓ,m \ m ∈ {0, . . . , n -1} 2 , ℓ = 0, . . . , L
is a redundant family of (L + 1)N vectors b Γ ℓ,m ∈ R N parameterized by the local direction flow Γ. Each grouplet vector b Γ ℓ,m is localized around the pixel m ∈ {0, . . . , n-1} 2 and follows the flow Γ on a length of 2 ℓ pixels. This family of vectors is thus efficient to compress a texture whose structures follow closely the direction of Γ. This is in particular the case for the image f 0 that has been used to extract the flow Γ as explained in section 2, but the tight frame B(Γ) can be used to analyze any image f of N pixels.

The forward grouplet transform [START_REF] Mallat | Geometrical grouplets[END_REF] corresponds to the computation of the decomposition of f onto the vectors of

B(Γ) = {b ℓ,m } ℓ,m ∀ ℓ, ∀ m, d ℓ (m) = f, b Γ ℓ,m . (2)
The grouplet vectors b Γ ℓ,m are defined implicitly through a decomposition algorithm that computes directly the coefficients d ℓ (m).

This decomposition is computed with a fast algorithm similar to the Haar wavelet transform. At a given scale ℓ, it processes the pixels m sequentially according to an ordering computed from the association field A ℓ , see [START_REF] Peyré | Texture synthesis with grouplets[END_REF]. It maintains a filtered residual f initialized to f . Each pixel m is linked to another pixel m = A ℓ (m). The grouplet detail coefficient d ℓ (m) is proportional to the difference f (m) -f ( m) and the value of the low pass residual f ( m) is updated to be proportional to the average ( f (m) + f ( m))/2. After processing the final scale ℓ = L -1, the low pass coefficients d L are set proportional to the low pass residual . The proportionality factors at each step of the transform are computed to ensure energy conservation even in the situation where several associations are converging from different points A ℓ (m 1 ) = A ℓ (m 2 ). The complexity of this fast grouplet transform is O(N (L + 1)) operations for an image of N pixels and L grouplet scales.

The backward grouplet transform retrieves an image f from the coefficients {d ℓ (m)} m,ℓ of the decomposition on the tight frame B(Γ). A stable inversion is obtained using the pseudo inverse that implements the reconstruction formula

f = 0 ℓ L τ ℓ m d ℓ (m)b Γ ℓ,m ,
where τ ℓ = 2 -ℓ for ℓ < L and τ L = 2 -L+1 . The fast backward grouplet transform performs this reconstruction and has a complexity of O(N (L+1)) operations, see [START_REF] Peyré | Texture synthesis with grouplets[END_REF].

Wavelet-Grouplet Transform

A grouplet basis vector b ℓ,m has an arbitrary length 2 ℓ but a fixed width of 1 pixel. This is problematic to represent textures with patterns of arbitrary width. Applying the grouplet transform over the wavelet coefficients of an image defines a waveletgrouplet transform. The corresponding wavelet-grouplet atoms have an arbitrary width 2 j that corresponds to the scale of the wavelet coefficients.

In this wavelet-grouplet setting, one first computes the decomposition of the image f on a set of wavelet vectors ∀ j = 0, . . . , J, ∀ p ∈ {0, . . . , n -1} 2 , f j (p) = f, ψ j,p , where ψ j,p is a wavelet vector scale 2 j and position p. A Laplacian pyramid tight frame is used in the numerical experiments, which has a redundancy of J + 1, see [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]. This multiscale transform has the advantage of being non-oriented, which leaves the processing of orientation to the grouplet transform alone. The reconstruction from the coefficients in such a tight frame reads

f = j κ j p f j (p)ψ j,p ,
where κ j = 2 -2j for j < J and κ J = 2 -2(J-1) . The set of coefficients f j is retransformed using the grouplet transform explained in section 2.2.

The projection of each f j on a grouplet frame B(Γ) generates the following set of coefficients

∀ j, ℓ, m, d j,ℓ (m) = f j , b Γ ℓ,m = f, b Γ j,ℓ,m
where b Γ j,ℓ,m is the wavelet-grouplet basis vector defined by

b Γ j,ℓ,m = p b Γ ℓ,m (p)ψ j,p ∈ R N .
A wavelet-grouplet atom b Γ j,ℓ,m is an elongated stroke of width ∼ 2 j and length ∼ 2 ℓ , that follows the geometrical flow Γ.

This set of vectors {b Γ j,ℓ,m } j,ℓ,m is a wavelet-grouplet tight frame B w (Γ) of R N composed of (J + 1)(L + 1)N vectors. The pseudo-inverse reconstruction from the wavelet-grouplet coefficients reads

f = j,ℓ τ j,ℓ m d j,ℓ (m)b Γ j,ℓ,m (3)
where τ j,ℓ = κ j τ ℓ .

The wavelet-grouplet forward transform corresponds to the computation of the coefficients {d j,ℓ (m)} j,ℓ,m of a given image f in the tight frame B w (Γ). This algorithm first computes the f j using the fast translation invariant wavelet transform, see [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF] and then applies the fast grouplet transform, to each of these f j . The fast backward wavelet-grouplet transform first applies the backward grouplet transform, and then the backward wavelet transform. Both forward and backward algorithms have a complexity of O((J + 1)(L + 1)N ) operations.

Texture Synthesis with Grouplets

Grouplet texture synthesis creates a new geometric texture f visually similar to a given input exemplar f 0 . This requires a statistical modeling of both the geometric layer and the wavelet-grouplet coefficients layer. The grouplet model for locally parallel textures is based on non parameteric statistics of the wavelet-grouplet coefficients layer. Similarly to the work of Heeger and Bergen [START_REF] Heeger | Pyramid-Based texture analysis/synthesis[END_REF], we retain only marginal statistics of the decomposition. In contrast to texture modeling over multiscale decompositions [START_REF] Heeger | Pyramid-Based texture analysis/synthesis[END_REF][START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF][START_REF] Zhu | Filters, random fields and maximum entropy (FRAME): Towards a unified theory for texture modeling[END_REF], our texture model also constrains the geometric layer of the texture.

Geometric and Coefficients Layers Learning

Grouplet layer learning. The geometric layer of the grouplet model is learned by computing the flow Γ 0 of the exemplar texture f 0 as detailed in Section 2.1.

Coefficients layer learning. The coefficients layer D 0 learned from f 0 is described through a the set of coefficients in the frame B w (Γ 0 )

∀ (j, ℓ), D j,ℓ (f 0 , Γ 0 ) = { f 0 , b Γ0 j,ℓ,m } m ,
together with the set of pixel values D sp (f 0 ) = {f 0 (x)} x . Each D j,ℓ (f 0 , Γ 0 ) is computed using the fast transform algorithm.

The layer D 0 is composed of several independent sets (one per scale (j, ℓ) and the pixel values) that empirical marginals of the statistical distribution of coefficients. They are used by the synthesis algorithm to enforce the synthesized texture to share the same empirical marginals as f 0 .

Coefficients Layer Synthesis

Texture synthesis with the grouplet model corresponds to synthesizing both a new geometric layer Γ and a new set of coefficients. This section assumes that the geometric layer Γ is known. The synthesized image is drawn at random from a texture ensemble parameterized by Γ. Section 4 describes how a dynamic layer can be computed using fluid dynamics evolutions.

Grouplet texture ensemble. Once both the geometric layer Γ 0 the coefficients layers D 0 have been learned, a grouplet texture ensemble T (Γ, D 0 ) is defined for any orientation field Γ. It is the set of textures having the same wavelet-grouplet and pixels marginals as f 0

T (Γ, D 0 ) = f ∈ R N \ ∀ (j, ℓ), D j,ℓ (f, Γ) = D j,ℓ (f 0 , Γ 0 ) ∩ f ∈ R N \ D sp (f ) = D sp (f 0 ) ,
This texture ensemble is thus parameterized by the geometric layer Γ of the new texture and by the coefficients layer D 0 of the exemplar. The geometry Γ that defines the model might be different from the learned geometry Γ 0 . This allows one to synthesize new textures with a different geometry.

Marginal constraint. Our texture model enforces synthesized coefficients to have the same coefficients layer as the one of the exemplar D 0 . Each marginal D j,ℓ,0 = D j,ℓ (f 0 , Γ 0 ) generates a constraint on the texture to synthesize. The synthesis requires to compute the orthogonal projection P(d j,ℓ , D j,ℓ,0 ) of wavelet-groulpet coefficients

d j,ℓ ∈ R N on d ∈ R N \ {d(m)} m = D j,ℓ,0 .
The same projection algorithm applies to the pixel values constraints generated by D sp (f 0 ). This projection corresponds to an histogram equalization of d j,ℓ with the values of D j,ℓ,0 , see [START_REF] Gonzales | Digital Image Processing[END_REF]. We denote by d r j,ℓ (i) the values of D j,ℓ,0 ordered by increasing values, and by d j,ℓ (r(i)) the ordered values of d j,ℓ , where r is a permutation of the indices. The projection is P(d j,ℓ , D j,ℓ,0 ) = dj,ℓ with dj,ℓ (r(i)) = d r j,ℓ (i).

If d j,ℓ and D j,ℓ,0 do not have the same number of coefficients, this formula requires interpolation.

Sampling the grouplet texture ensemble. As explained in the FRAME model [START_REF] Zhu | Filters, random fields and maximum entropy (FRAME): Towards a unified theory for texture modeling[END_REF], the uniform distribution on T (Γ, D 0 ) has maximal entropy, thus leading to the least synthesis bias. Performing an un-biased synthesis is difficult, and Portilla and Simoncelli [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF] replace this uniform sampling by a sampling with high entropy. This is achieved by iteratively projecting a white noise image on the set of constraints that define the texture ensemble. The grouplet coefficients synthesis algorithm detailed in Table 1 uses a similar iterative projection strategy. The initial image f (0) is a random Gaussian white noise image, and the synthesis algorithm converges to a synthesized image S(f (0) ; Γ, D 0 ) ∈ T (Γ, D 0 ). The procedure iteratively decomposes the image into the wavelet-grouplet frame B w (Γ), forces the wavelet-grouplet marginal D j,ℓ (f, Γ) to match the coefficients layer D j,ℓ,0 of the exemplar, reconstructs an image from the grouplet coefficients and then enforces consistency of the pixel marginal D sp (f ) with the one of the exemplar.

Figure 2 shows the iterations of the synthesis with a user defined geometric layer Γ. The iterations progressively filter the noise along the geometric flow to create texture patterns. The exemplar image f 0 is twice smaller than the synthesized texture f . Input: orientation flow Γ, initial image f (0) , coefficients layer D 0 = {D j,ℓ,0 } j,ℓ ∪ D sp 0 . Output: synthesized texture f .

(1) Initialization: set i ← 0.

(2) Forward transform: compute the wavelet-grouplet coefficients {d

(i) j,ℓ (m)} j,ℓ,m of f (i) in B w (Γ).
(3) Grouplet constraints: for each (j, ℓ), perform the histogram equalization (4): d(i) j,ℓ ← P(d

(i) j,ℓ , D j,ℓ,0 ). (4) Backward transform: compute f (i+1) from the coefficients {d (i) j,ℓ (m)} j,ℓ,m .
(5) Projection on gray-level constraint: equalize the pixel histogram: f (i+1) ← P( f (i+1) , D sp 0 ). (6) Stop: while not converged, set i ← i + 1 and go back to 3.

Table 1: Synthesis algorithm to compute S(f (0) ; Γ, D 0 ).

Dynamic Texture Synthesis

A dynamic fluid texture { ft } t 0 is created from a single exemplar image f 0 by evolving in time the direction field Γ t a of the geometric layer according to the incompressible Navier Stokes equation

∂Γ t ∂t = P inc (-(Γ t • ∇)Γ t + ν∆Γ t ) ,
with Neumann or periodic boundary conditions, where ν is the viscosity of the fluid and ∆ is the Laplacian operator. The projector ṽ = P inc (v) on the set of divergence free vector fields, computed by solving the Poisson equation for a scalar potential V ṽ = v -∇V where ∆V = ∇ • v.

The synthesized texture ft for t > 0 is obtained by solving an advection equation projected on the grouplet texture model T (Γ t , D 0 ) defined in Section (3.2) by the current geometry Γ t and the coefficients layer D 0 of f 0

∂ ft ∂t = S -(Γ t • ∇) ft + ν∆ ft ; Γ t /||Γ t ||, D 0 , (6)
where ν is the viscosity of the texture, that might be different from ν. The projection S(• ; Γ, D 0 ) on the grouplet model is computed with a few steps of the iterative algorithm detailed in Table 1.

The solution of both ( 5) and ( 6) are computed numerically using the implicit solver of [START_REF] Stam | Stable fluids[END_REF]. It makes use of the following flow-warping operator, that is defined for scalar function and vector fields W Γ (f ) = f where f (x) = f (x -Γ(x)). [START_REF] Efros | Texture synthesis by non-parametric sampling[END_REF] Table 2 details the steps of the fluid texture synthesis algorithm 1 . Figure 3 shows two examples of dynamic texture synthesis.

Input: initial texture f 0 , time step size τ > 0.

Output: synthesized texture fi for times 0 τ i T .

(1) Initialization: compute Γ 0 and D 0 from f 0 , set i ← 0.

(2) Advect: compute Γ

= W τ Γi (Γ i ) and f

= W τ Γi ( fi ).

(3) Diffuse: compute

Γ (2) i = Γ (1)
i + τ ∆Γ Table 2: Grouplet fluid texture synthesis.

Conclusion

This paper proposed a new algorithm to perform fluid texture synthesis. The analysis of the exemplar texture is carried over by computing a geometric flow that drives the wavelet-grouplet transform. This defines a texture model using statistics of the wavelet-grouplet coefficients of the exemplar. A dynamic texture is obtained by 1 The turbulence image is obtained by the resolution of fluid dynamics equations computed by M. Farge et al. [START_REF] Farge | Wavelets and turbulence[END_REF]. The second texture is obtained from [START_REF] Liang | Real-time texture synthesis by patch-based sampling[END_REF]. synthesizing the geometry using fluid dynamics and by evolving in time the texture inside the texture model. This opens many questions regarding the dependancy of geometry and texture in natural images, which is only partially addressed by our model.

Figure 1 .

 1 Figure 1. Examples of locally parallel textures.

Figure 2 .

 2 Iterations of the synthesis algorithm 1.

( 4 )( 5 )

 45 Project: compute Γ i+1 = P inc (Γ i+1 /||Γ i+1 ||, D 0 ) Stop: while τ i < T , go back to 2.

Figure 3 .

 3 Example of dynamic texture synthesis. The exemplar images used are the first two textures shown on Figure 1.

CNRS and Ceremade, Université Paris-Dauphine, 75775 Paris Cedex 16 France. c PUO -Septembre

Acknowledgements. I would like to thank Stéphane Mallat and Eero Simoncelli for fruitful discussions, and Marie Farge for providing the turbulence image.