Ivan Enderlin
email: ivan.enderlin@femto-st.fr

Alain Giorgetti
email: alain.giorgetti@femto-st.fr

Fabrice Bouquet
email: fabrice.bouquet@femto-st.fr

A Constraint Solver for PHP Arrays

Keywords: array, constraint, solver, php, realistic domain, praspel

In previous works, we have proposed Praspel, a framework for contract-based testing in PHP. Among others, it includes a specification language and a unit test generator which automatically generates test data from formal preconditions. The generator sometimes rejects data, when they do not satisfy parts of the preconditions. In many cases, generation with rejection is not efficient enough. Thus we investigate practical contexts where more efficient generation algorithms can be designed and we extend Praspel with their implementation.

After strings, that we have already considered, the most frequent data type in PHP is arrays. They cover most of the needs for collections, because they can store key-value pairs of any kind, they do not have a specific length or depth, and they are efficiently implemented. In this paper, we report on a study to know what are the most popular constraints on PHP arrays. Then we formalize these constraints and we present an implementation in PHP of a constraint solver for these constraints. In this context, the constraint-based approach removes all the rejections.

I. INTRODUCTION Contract-based testing

 [START_REF] Aichernig | Contract-based testing[END_REF]is a promising approach to increase software safety and security. It is based on the notion of Design by Contract, introduced by Meyer [START_REF] Meyer | Applying "design by contract[END_REF] with the Eiffel language [START_REF]Eiffel: programming for reusability and extendibility[END_REF]. A contract is a formal specification associated to the code of a program. It mainly consists of invariants, preand postconditions. Invariants describe properties that should hold at each step of the execution. Pre-and postconditions respectively represent conditions that have to hold for a method to be invoked, and conditions that have to hold after execution of the method.

Various contract languages extend programming languages with contracts: JML extends Java [START_REF] Leavens | JML: A notation for detailed design[END_REF], ACSL extends C [START_REF] Baudin | ACSL: ANSI C Specification Language[END_REF], Spec# extends C# [START_REF] Barnett | The Spec# Programming System: An Overview[END_REF], etc. Among the numerous advantages of contracts, formal properties can be exploited for (unit) testing. Indeed, the information contained in invariants and preconditions can be used to generate test data. In addition, these assertions can be checked at run time and thus provide a (partial) test oracle for free.

In previous works, we have introduced Praspel [START_REF] Enderlin | Praspel: A specification language for contract-based testing in php[END_REF], a toolsupported specification language for contract-based testing in PHP [START_REF] Zaïdi | The PHP website[END_REF]. Praspel extends contracts with the notion of realistic domain, which makes it possible to assign a domain of values to data (class attributes or method parameters). A library of predefined realistic domains is already available along with a test environment [START_REF] Enderlin | Hoa project, a set of PHP libraries[END_REF].

Validating PHP web applications often involves to manipulate strings and arrays. Strings have been already ad-dressed [START_REF] Enderlin | Grammar-Based Testing Using Realistic Domains in PHP[END_REF] by introducing grammar-based testing in Praspel. This paper focus on PHP arrays, which are always associative arrays. They cover most of the needs for collections, because they can store key-value pairs of any kind (they can be homogeneous and heterogeneous), they do not have a specific length or depth, and they are efficiently implemented. When generating arrays for testing, the main difficulty is to satisfy their properties, which are formalized by predicates in their specification. We address this issue in the context of PHP and Praspel, by designing a specific constraint solver.

Our contributions are twofold. First, after studying the most popular conditions on PHP arrays, we propose a syntax to specify them in Praspel and a semantics of these array specifications. Second, we propose a new constraint solver in PHP for PHP, to generate arrays satisfying these conditions. These two contributions are embedded in realistic domains and can be used in Praspel annotations of a PHP program.

The paper is organized as follows. Section II briefly introduces the notion of realistic domain and its implementation in Praspel for PHP. Section III presents a language of array conditions in Praspel, inspired by a study to know what are the most popular conditions expressed on arrays in popular PHP code. Section IV presents the constraint solver. Then, Section V reports on an experimentation validating our approach and tool, and showing its usefulness and efficiency in practice. Related works are presented in Section VI. Finally, Section VII concludes and presents future works.

II. REALISTIC DOMAINS AND PRASPEL

This section is a reminder of [START_REF] Enderlin | Praspel: A specification language for contract-based testing in php[END_REF]. It shortly presents the notion of realistic domain and its application to PHP programs. Realistic domains are designed for test generation purposes. They specify which values can be assigned to a data in a given program. They are well-suited to PHP, since this language is dynamically typed (i.e. no types are syntactically assigned to data) and realistic domains thus introduce a specification of data types mandatory for test data generation. We first introduce general features of realistic domains. Then we present their implementation in PHP and the Praspel framework.

A. Features of Realistic Domains

Realistic domains refine usual datatypes (integer, string, array, etc.). They are intended to specify data domains relevant for specific application contexts. For example, email addresses constitute a realistic domain: many applications identify a user by her email address, so we need to validate and generate such data, and an email address is more than a string: it obeys to a specific syntax that makes it not obvious to generate.

The first feature of a realistic domain (named its predicability) is to carry a characteristic predicate of its values, used to check if a value belongs to the possible set of values described by the realistic domain. The second feature of a realistic domain (named its samplability) is to propose a value generator, called the sampler, that generates values in the realistic domain. For instance, a sampler for the realistic domain of email addresses can generate strings representing syntactically correct email addresses defined by a regular expression.

B. Realistic Domains in PHP

In PHP, we have implemented realistic domains as classes providing at least two methods, corresponding to the two features of realistic domains. The first method is named predicate($q) and takes as input a value $q: it returns a boolean indicating the membership of the value to the realistic domain. The second method is named sample() and generates values that belong to the realistic domain.

Our implementation of realistic domains in PHP exploits the PHP object programming paradigm and takes benefit from the following three principles:

1) Inheritance: PHP realistic domains can inherit from each other. A realistic domain child inherits the two features of its parent, namely predicability and samplability, and is able to redefine them. Consequently, we say that all the realistic domains constitute a hierarchical universe.

2) Interfaces: All data are represented by realistic domains. Some of them implement interfaces which characterize them. Useful interfaces for our current concerns are:

• Constant: represents an immutable realistic domain with only one value, such as 42, true, etc. • Interval: represents an interval by its lower and upper bounds, that can be dynamically reduced, • Nonconvex: allows to discredit values, i.e. to specify that a value no longer belongs to a realistic domain and should therefore not be generated, • Finite: allows to count the number of values, • Enumerable: allows to iterate over all the values. Thus, a realistic domain that implements the Interval and Nonconvex interfaces is an interval with "holes". Counting and exhaustive generation of a finite realistic domain take discredited values into account.

3) Parameterization: Realistic domains may have parameters, like a function does. Data given to a realistic domain are called arguments. This feature helps generating structured data such as arrays, objects, graphs, automata, etc.

Example 1 (Realistic domains with arguments). The realistic domain string(0x61, 0x7a, boundinteger [START_REF] Leavens | JML: A notation for detailed design[END_REF][START_REF] Cheon | A Simple and Practical Approach to Unit Testing: The JML and JUnit Way[END_REF]) admits as arguments two integers (that represent two Unicode code-points) and a domain of integers to specify its possible length. The realistic domain boundinteger(X, Y ) contains all the integers between X and Y . We can also write X..Y as syntactic sugar for this domain. The realistic domain string(X, Y , L) is intended to contain all the strings of length (in the domain) L built of characters from X to Y code-points.

When describing parameters, because all data are realistic domains, we can multi-type hint the parameters, i.e. filter with multiple realistic domain names. For instance, the first two parameters of the realistic domain string() are described as Constinteger | Conststring, so can be a constant integer or a constant string. Then, the syntax string('a', 'z', 4..12) is strictly equivalent to the previous one. If we write string(true, 'z', 4..12), an error will be thrown. The realistic domain itself handles the cast.

C. Praspel and Contract-Based Testing in PHP

Praspel means PHP Realistic Annotation and SPEcification Language. It is a language and a framework for contract-based testing in PHP, based on realistic domains.

Praspel annotations are written inside comments in the source code. Invariants document classes and pre-and postconditions document methods. The general form of Praspel annotations is shown in Figure 1. For lack of space, named behaviors and specification of exceptions are not presented. In this specification, I 1 , . . . , I h are invariant clauses, assumed to be satisfied at the beginning and at the end of each method invocation. Formulas R 1 , . . . , R n and A 1 , . . . , A k are precondition clauses, that have to hold at the invocation of the foo method. Formulas E 1 , . . . , E m are postconditions that have to be established when method foo terminates without throwing an exception. In postconditions, Praspel provides the two additional constructs \result and \old(e), which respectively designate the value returned by the method, and the value of expression e at the pre-state of the method invocation.

PHP does not provide a type system, but Praspel contracts make it possible to give typing informations, assigning realistic domains to data (class attributes or method parameters). The construction i: t 1 (. . . ) or . . . or t n (. . . ) associates at least one realistic domain, among t 1 (. . . ), . . . , t n (. . . ), to the identifier i. Every identifier holds a realistic domain disjunction, thanks to the or keyword. Identifiers can also be passed into arguments of realistic domains.

Example 2 (Identifiers in Praspel). The declarations:

length: 4..12 str : string('a', 'z', length)

show a dependency between two identifiers: The third argument of the realistic domain of the second identifier str is the first identifier length.

Praspel provides a set of more than 30 predefined realistic domains, called the standard library. Some of them correspond to scalar types (integer, float, boolean. . . ), other ones to classes and arrays (detailed in Section III).

Contractual assertions are made of realistic domain declarations, possibly completed with additional predicates, expressed in PHP using the \pred construct.

III. ARRAYS IN PRASPEL

In PHP, an array is always an associative array (or map, or dictionary), i.e. a collection of key-value pairs, where each key appears at most once. Keys can be null, booleans (casted into integers), integers, floats (reduced to their integer parts) or strings. Values can be of many kinds. An array can be homogeneous or heterogeneous. In an homogeneous array all the keys have the same type, and all the values too. In an heterogeneous array keys may have distinct types, and/or values may have distinct types. Keys can be auto-incremented, by adding 1 to the last integer index starting by 0. The length (or size) of an array is its number of elements. An array has no predefined length, but its length (stored internally by the PHP engine) can be retrieved thanks to the PHP function count(). An array has also no predefined depth, i.e. it can contain arbitrary arrays.

A. Array Description

In Praspel, array(D,L) denotes the realistic domain of arrays whose domains and codomains are described by D and whose length is in the disjunction L of realistic domains of non-negative integers. D is a comma-separated list, between [ and ], of array descriptions of the form from K to V , where K and V are realistic domain disjunctions, respectively for keys and values. When the from keyword is missing, it is transformed into a realistic domain representing an autoincremented integer starting at 0 with a step of 1. The identifier a1 is declared as a homogeneous array of booleans with a length between 7 and 42. This length is a realistic domain that implements the Interval interface. The identifier a2 is declared as a homogeneous array of length 7, whose keys are integers between 0 and 5 or simply 10, and whose values are integers. Its length is a realistic domain that implements the Constant interface, the domain for its keys is the disjunction of two realistic domains (0..5 and 10..23) which implement the Interval interface, etc. The identifiers a3 and a4 are declared as heterogeneous arrays. Both arrays can contain the pairs (5, true) and (15, 4.2), but a4 can contain the pair (5, 4.2), whereas a3 can not contain it.

Example 3 (Homogeneous and heterogeneous arrays). This syntax is illustrated by the following array declarations:

Actually, we introduce a normal form that removes disjunctions in array descriptions (in from . . . to . . . constructs). An array description is in normal form when it can not be reduced by the rewriting rule (from

F 1 or F 2 to T 1 or T 2 → from F 1 to T 1 , from F 1 to T 2 , from F 2 to T 1 , from F 2 to T 2 ).
Example 4 (Array description in normal form). The following declaration of a4 is in normal form:

B. Collecting Information

When generating test data, the Praspel testing tool calls the sample() method on a realistic domain and then checks whether the predicates declared with \pred() hold. For an array with conditions, a random generator may produce a lot of rejected data before getting a valid one. The idea is to determine the most popular conditions on arrays expressed in PHP (usually written in the \pred() construct), to allow them inside Praspel, and to use a solver to satisfy these conditions.

To achieve this, we have selected 61 PHP projects, from Github and SourceForge, for their popularity, impact on the industry and complexity. All these projects represent 28 066 files and 5 220 547 lines of code. In this code we count the number of occurrences of each array function available in the PHP standard library. It appeared that the three most used functions are: count(), array_key_exists() and in_array(). The count() function counts the number of values in an array, the array_key_exists() function checks whether a key is present in an array (independently of its associated value, e.g. it returns true even if the value is null), and finally, the in_array() function checks whether a value is present in an array. All these functions work on one array at a time. This study suggests that we could consider these side-effect-free Boolean functions as the most frequent conditions on arrays.

C. Array Conditions

We extend the syntax of the array declaration a : array(D,L) in Praspel with the following conditions on arrays.

A pair condition is of the form a[K]: V where K and V are realistic domain disjunctions. The condition means that the pairs constituted of all the keys in K and at least one value in V are present in the array a. K only accepts realistic domains that implements the Constant, Interval and Enumerable interfaces. This is equivalent to use the array_key_exists() and in_array() functions combined.

If we would like to express a constraint only on K, we can use the symbol _. The condition a[K]: _ means that all the keys from K must be present in the array a. It is equivalent to use only the array_key_exists function with all values in K in conjunction.

The condition a[_]: V means that all the values in V must be present in the array a. It is equivalent to use the in_array function with all the values in V in conjunction.

Instead of the : symbol, we can use the symbol !: to express a negation. The condition a[K]!: V means that all the keys in K have a value in the array a and that this value is not in V . It works similarly with the symbol _. For example, the condition a[K]!: _ means that no key in K appears in a.

Keys of an array are always unique, but not its values. We can express a unicity condition on values by writing the condition a is unique. In this case, we cannot have the same value twice in the array a.

Example 5 (Array Conditions).

To illustrate all kinds of conditions, we will use the following example that uses a: length: 0.. [START_REF] Baudin | ACSL: ANSI C Specification Language[END_REF]

or 10 a : array([to string('a', 'e', 1)], length) a[0]: 'b' or 'd' a is unique

IV. CONSTRAINT SOLVER

Given a conjunction of array conditions on an array a, we propose to invoke a constraint solver to construct an array satisfying all these conditions. This section explains how array conditions are transformed into constraints for the solver. One of these conditions is assumed to be an array declaration of the form a: array(D,L), where D is assumed to be in normal form, without loss of generality. In other words, D is a list of p constructs from F i to T i with 1 ≤ i ≤ p. We also assume that L is L 1 or . . . or L m (with m ≥ 1), where L 1 , . . . , L m are realistic domains that inherit from the Integer realistic domain and that are non-negatives.

In Example 5, p = 1, m = 2, L 1 = [0..5] and L 2 = {10}. In the array description, no domain is declared. In this case the default realistic domain integerpp(0,1) is used (an autoincremented integer: 0, 1, 2, 3, etc). In this example F 1 = integerpp(0,1) and T 1 = string('a','e',1).

Without risk of confusion, a domain disjunction D 1 or . . . or D n will often be identified with the set

D 1 ∪ • • • ∪ D n .

A. Variables

The constraint variables are: (i) the array size -noted Swhich is a non-negative integer, (ii) the sets X and Y , which respectively are the array domain (set of keys) and codomain (set of values), (iii) the array content noted H, which is a total function from X to Y , since keys are unique in a PHP array, (iv) the realistic domains 1 X 1 , . . . , X p (resp. Y 1 , . . . , 1 In fact, the solver will not handle realistic domains, but only sets.

Y p ), which are subsets of the realistic domains F 1 , . . . , F p (resp. T 1 , . . . , T p ) compatible with all the array conditions.

We are essentially interested in finding the content of X and the values of the function H, i.e. the content of H considered as a hashtable, possibly also the values of the X i s and Y i s for checking purposes. The other variables are introduced only to simplify the expression of constraints.

When x ∈ X holds, H(x) = y means that the key-value pair (x, y) is in the array. We extend H to subsets of X by the function Ĥ defined by Ĥ(E) = {H(x) s.t. x ∈ E} for any subset E of X.

B. Cardinality Constraints

Let card(E) denote the cardinality of the finite set E. The constraints card(X) = S and S ≥ 0 say that the array size is its number of keys and is non-negative.

By default, there is no unicity constraint on the codomain, so we only have the constraint card(Y ) ≤ card(X) However, in presence of the array condition a is unique, this constraint becomes card(Y ) = card(X).

C. Constraints on the Array Size

When propagating constraints, the solver may refine the domains L 1 , . . . , L m of possible values for the array size S for i in {1, . . . , m}, and the array size S should be in one of these sets, i.e. the constraint to define domain of S is:

S ∈ L 1 ∪ • • • ∪ L m In Example 5, we have L 1 ⊆ [0..5] and L 2 ⊆ {10}. The size S is constrained by S ∈ L 1 ∪ L 2 .

D. Constraints on Domains and Codomains

The domain X and codomain Y of H are related by the constraint Y = Ĥ(X).

We expect that the constraint solver proposes us the array domain X (resp. codomain Y ) as a disjunction X 1 or . . . or X p (resp. Y 1 or . . . or Y p ) of realistic domains compatible with all the array conditions. We should have the equalities X = 1≤i≤p X i and Y = 1≤i≤p Y i , and the inclusions:

X i ⊆ F i and Y i ⊆ T i for i in {1, . . . , p}. The pair (X i , Y i ) with 1 ≤ i ≤ p should also satisfy the constraint Ĥ(X i ) = Y i meaning that Y i is the codomain of the restriction of H to X i (⊆ X).

E. Constraints on Pairs

For each array condition a[K]: V where K and V are domain disjunctions, we introduce the constraints: K ⊆ X and Ĥ(K) ⊆ V . A negated pair condition a[K]!: V is translated into the constraints: K ⊆ X and Ĥ(K) ∩ V = ∅.

F. Constraints on Keys or Values

The condition a[K]: _ is translated into the constraint K ⊆ X, and its negation a

[K]!: _ into the constraint K ∩ X = ∅. The condition a[_]: V is translated into the constraint V ⊆ Y , and its negation a[_]!: V into the constraint V ∩ Y = ∅.

G. Propagation and consistency

Propagation of constraints uses an AC3 algorithm [START_REF] Macworth | Consistency in network of relations[END_REF] implemented in PHP. So, we use five kinds of domains associated to five kinds of realistic domains: Constant, Interval, Nonconvex, Finite and Enumerable (see Section II-B2). For each kind of domain, we have implemented a revise method to allow the domain reduction. So, the consistency is also checking that there is no empty domain for the four variables S, H, X and Y but not for X i and Y i . The goal is to detect inconsistencies as soon as possible.

H. Labelling

The labelling is the process of finding a value for each variable. In order to make the solver converge quickly to a solution, we use a heuristic that consists in chosing a value for the variable S at first. This helps to unfold the ∀ and ∃ quantifiers (because F i and T i are enumerables, we manipulate finite sets). Then, the solver tries to compute the sets X i and Y i . We use a random generator to generate a value in a realistic domain, to select a realistic domain in a disjunction, etc. The generated value is then propagated. If an inconsistency is detected, we add a new constraint to discredit the value, and then generate another one. For instance, if S = 5 leads to an inconsistency, we add the constraint S = 5. The added constraint is removed during the backtracking step.

When all variables are labelled, i.e. each one has a valid value, the solver returns the solution.

V. EXPERIMENTATION

In all that follows, the word system stands for the expression "conjunction of array conditions". This section presents an experimentation evaluating the solver efficiency, i.e. its capability to avoid or reduce rejection when generating data from systems. We measure the number of backtracks in the solver, the time to generate data from satisfiable systems of array conditions, and how many unsatisfiable systems are detected. The experimentation is composed of three steps: system generation, then data generation (i.e. system solving) and finally a measuring step. We generate systems on arrays containing strings and integers, and of length 5 to 20.

The Praspel language (and in particular its sublanguage of array conditions) is described by a grammar. In order to generate systems, we re-use a previous work in grammar-based testing [START_REF] Enderlin | Grammar-Based Testing Using Realistic Domains in PHP[END_REF] proposing three algorithms generating data from grammars: a uniform random generator, a bounded exhaustive test generator, and a rule-coverage-based test generator. Since the grammar of array conditions is small, the last generator does not generate a sufficiently wide collection of systems. Bounded exhaustive testing is more costly than random testing, but it is more precise and well adapted to small grammars. We retain the bounded exhaustive test generator: for increasing values of n it enumerates all the systems composed of a sequence of n tokens: In this generator, a single token value of each variable token (i.e. interval bounds and particular array lengths, keys and values) is generated at random.

With n = 3, we can generate conditions of the form arr is unique. With n = 6, we can generate constraints of the form arr[0]: 0. With n = 8, we can generate constraints of the form arr[0]: 0 or 1. With n = 11, we can generate for instance the system in Example 5.

The second step calls the solver with each produced system to generate an array satisfying it. Every generated array is evaluated by the predicate associated to the system of array conditions, to check the solver soundness. During the data generation step, some marks are placed to count only the array generation time, without counting the compiling time. We also measure the number of backtracks in the solver and the number of rejected systems. When a system is rejected, its backtracks are not counted.

Since our grammar-based testing algorithms use an isotropic random generation to generate token values, these values can differ from one system to another. It may lead to different generation times and numbers of backtracks. To avoid peaks in the results, we report the average of 100 generations of systems sharing the same pattern.

Table I shows our experimentation results. In the first column, n is the length of the generated sequences of tokens representing systems of array conditions. Column 2 gives the number of distinct system patterns with this length. Columns 3, 5 and 6 respectively give average numbers of backtracks and rejected systems for 100 runs, and an average generation time. Column 4 gives the rate of backtracks per system pattern.

For n ≤ 20, we observe that no system is rejected, which is a great improvement by comparison with random generation. All generated data satisfy their specification. The solver successfully and quickly generates data with a low number of backtracks. For n = 20, which represents approximatively 3 constraints with disjunctions, the exhaustive generation algorithm generates 492 distinct systems of array conditions and the execution produces only 114 backtracks, so approximatically 1 backtrack for 4 systems. This is a good result. Nevertheless, we were not able to characterize the number of backtracks with the number of constraints. Finally, this experimentation allowed us to find a bug in our solver. This bug always led to a rejection when only a certain constraint was analyzed. Thus, it shows that the validation process also helps finding bugs.

VI. RELATED WORKS

Various works consider Design-by-Contract for unit test generation [START_REF] Cheon | A Simple and Practical Approach to Unit Testing: The JML and JUnit Way[END_REF], [START_REF] Madsen | Unit Testing using Design by Contract and Equivalence Partitions[END_REF]. Our specification language is based on JML [START_REF] Leavens | JML: A notation for detailed design[END_REF] and ACSL [START_REF] Baudin | ACSL: ANSI C Specification Language[END_REF]. Thanks to an expressive specification language, Praspel performs general runtime assertion checks. Realistic domains present some similarities with Eiffel's types [START_REF]Eiffel: programming for reusability and extendibility[END_REF], especially regarding inheritance between realistic domains. Nevertheless, the two properties of predicability and samplability displayed by realistic domains do not exist in Eiffel. Moreover, Praspel adds clauses that Eiffel contracts do not support, as @throwable and @behavior, which are inspired from JML.

Euclide [START_REF] Gotlieb | Euclide: A Constraint-Based Testing Framework for Critical C Programs[END_REF] is a constraint-based testing tool that could take as input additional safety properties defined in ACSL [START_REF]Tcas software verification using constraint programming[END_REF]. Our approach differs by handling conditions directly in the contract. All constructions present in Praspel are well-handled for both aspects: validation and generation. The CLP framework INKA [START_REF] Gotlieb | A CLP Framework for Computing Structural Test Data[END_REF] helps computing structural test data from a C program. It transforms the problem of automatic test data generation into a CLP problem over finite domains. In the same way, FDCC [START_REF] Bardin | FDCC: A Combined Approach for Solving Constraints over Finite Domains and Arrays[END_REF] is a combined approach for solving constraints over finite domains and arrays. The tricky part of FDCC lies in a bi-directional communication mechanism between two solvers. Our constraints are more specific but we use only one solver.

VII. CONCLUSION AND FUTURE WORKS

We have presented in this paper an extension of the Praspel language to specify usual conditions on PHP arrays. We have expressed its semantics by constraints. We have designed and implemented in PHP a constraint solver to generate test data from these constraints. It uses a random generator to ensure a diversity of generated solutions. This solver is integrated in realistic domains and can be used within the Praspel framework. A first validation shows cases where rejection has been totally removed. It also shows that the solver dramatically increases the generation speed.

In a near future, we plan to lead a more complete experimentation. Then, we plan to formalize more constraints and extend our solver. We also plan to transform constraints into the formalism proposed by MiniZinc [START_REF] Nethercote | Minizinc: Towards a standard cp modelling language[END_REF] in order to compare our solver to other ones regarding performances and capabilities to find solutions. As a string is an array of characters, we would like to apply the same process on strings with the help of existing and promising solving techniques [START_REF] Ganesh | Hampi: A string solver for testing, analysis and vulnerability detection[END_REF].

  a1: array([to boolean()], 7..42) a2: array([from 0..5 or 10 to integer()], 7) a3: array([from 0..10 to boolean(), from 20..30 to float()], 7) a4: array([from 0..10 or 20..30 to boolean() or float()], 7)

  a4: array([from 0..10 to boolean(), from 0..10 to float(), from 20..30 to boolean(), from 20..30 to float()], 7)

  For the condition a[0]: 'b' or 'd' in Example 5, we have K = {0} and V = {'b', 'd'}. The constraints are {0} ⊆ X and Ĥ({0}) ⊆ {'b', 'd'}.

  class C { / ** @invariant I 1 and . . . and I h * / / ** @requires R 1 and . . . and Rn; * @ensures E 1 and . . . and E j ; * / function foo ( $x 1 . . . ) { body } } Figure 1. Syntax of contracts in Praspel Test generation in Praspel is decomposed into two steps. First, a test generator computes test data from contracts. Second, a dedicated test execution framework runs the test cases (i.e. invokes the methods with the computed test data, and checks the assertions at run time) so as to establish the test verdict.


Patrick Héas
email: patrick.heas@inria.fr

Frédéric Lavancier

Souleymane Kadri-Harouna

P H Éas

P D Érian

Consider instead the following article: "Self-similar prior and wavelet bases for hidden incompressible turbulent motion", by

Keywords: Bayesian modeling, fractional Brownian motion, divergence-free wavelets, fractional Laplacian, wavelet connection coefficients, fast recursive filtering, optic-flow, image restoration AMS subject classifications. 60G18, 60G22, 60H05, 62F15, 65T50, 65T60

This paper presents a method for regularization of inverse problems. The vectorial bi-dimensional unknown is assumed to be the realization of an isotropic divergence-free fractional Brownian Motion (fBm). The method is based on fractional Laplacian and divergence-free wavelet bases. The main advantage of these bases is to enable an easy formalization in a Bayesian framework of fBm priors, by simply sampling wavelet coefficients according to Gaussian white noise. Fractional Laplacians and the divergence-free projector can naturally be implemented in the Fourier domain. An interesting alternative is to remain in the spatial domain. This is achieved by the analytical computation of the connection coefficients of divergence-free fractional Laplacian wavelets, which enables to easily rotate this simple prior in any sufficiently "regular" wavelet basis. Taking advantage of the tensorial structure of a separable fractional wavelet basis approximation, isotropic regularization is then computed in the spatial domain by low-dimensional matrix products. The method is successfully applied to fractal image restoration and turbulent optic-flow estimation.

Introduction. Fractional Brownian motion (fBm) introduced by Mandelbrot and Van Ness

 [START_REF] Mandelbrot | Fractional brownian motions, fractional noises and applications[END_REF]constitute popular models for the study of self-similar stochastic phenomena in a broad range of application domains such as statistical turbulence, network traffic analysis, financial markets or image processing. The practical interest of fBm models relies in their remarkable properties: Gaussianity, increments stationarity and statistical self-similarity. The latter property goes hand-in-hand with the empirical observation of inverse power-law power spectra in a wide range of applications, although the latter statement has to be made with caution as classical notion of power spectrum are not applicable for non-stationary signals [START_REF] Flandrin | On the spectrum of fractional brownian motions[END_REF]. Extensions of fBm have been proposed in the context of multidimensional signal analysis. In particular, there has been a special focus on the isotropic case [START_REF] Reed | Spectral representation of fractional brownian motion in n dimensions and its properties[END_REF], since isotropy constitutes a simple, yet relevant, assumption in numerous domain, including statistical turbulence [START_REF] Monin | Statistical Fluid Mechanics: Mechanics of Turbulence[END_REF]. Finally, vectorial fBm models first proposed in [START_REF] Monin | Statistical Fluid Mechanics: Mechanics of Turbulence[END_REF] have been recently studied systematically in [START_REF] Amblard | Basic properties of the multivariate fractional Brownian motion[END_REF][START_REF] Tafti | Fractional Brownian vector fields[END_REF].

Interesting fBm decomposition can be designed by projection of fBm on wavelet bases. Indeed, non-stationarity and self-similarity of fBm necessitate some appropriate space-scale analysis brought by wavelet multiresolution analysis [START_REF] Flandrin | Wavelet analysis and synthesis of fractional brownian motion[END_REF]. In particular, bi-orthogonal fractional wavelets constitute whitening filters for these signals, which are correlated both in space and scale [START_REF] Elliott | A fourier-wavelet monte carlo method for fractal random fields[END_REF][START_REF] Meyer | Wavelets, generalized white noise and fractional integration: the synthesis of fractional brownian motion[END_REF][START_REF] Tafti | Invariances, laplacian-like wavelet bases, and the whitening of fractal processes[END_REF][START_REF] Tafti | Fractional Brownian vector fields[END_REF]. Moreover, in the context of turbulence studies, bi-orthogonal wavelet bases can also be designed to represent the space of divergence-free functions, which is a constraint of the incompressible Navier-Stokes equations [START_REF] Lemarié-Rieusset | Analyses multirésolutions non orthogonales, commutation entre projecteurs et dérivation et ondelettes vecteurs à divergence nulle[END_REF]. Another advantageous property is the existence of a fast implementation of decomposition-reconstruction algorithms for those wavelets using recursive filter banks [START_REF] Abry | The wavelet-based synthesis for fractional brownian motion proposed by f. sellan and y. meyer: Remarks and fast implementation[END_REF][START_REF] Deriaz | Divergence-free and Curl-free wavelets in 2D and 3D, application to turbulence[END_REF].

FBm models have intensively been studied in the context of signal synthesis or analysis. In particular, many works deal with the synthesis of fBm (see e.g., [START_REF] Abry | The wavelet-based synthesis for fractional brownian motion proposed by f. sellan and y. meyer: Remarks and fast implementation[END_REF][START_REF] Amblard | Basic properties of the multivariate fractional Brownian motion[END_REF][START_REF] Flandrin | Wavelet analysis and synthesis of fractional brownian motion[END_REF][START_REF] Elliott | A fourier-wavelet monte carlo method for fractal random fields[END_REF][START_REF] Meyer | Wavelets, generalized white noise and fractional integration: the synthesis of fractional brownian motion[END_REF][START_REF] Tafti | Fractional Brownian vector fields[END_REF]) or with the characterization of Hurst exponent from experimental data (see e.g., [START_REF] Flandrin | Wavelet analysis and synthesis of fractional brownian motion[END_REF][START_REF] Blu | Self-similarity: Part ii optimal estimation of fractal processes[END_REF][START_REF] Tafti | Invariances, laplacian-like wavelet bases, and the whitening of fractal processes[END_REF]). However, suprinsingly, very little studies concern Bayesian analysis and fBm models, i.e., employing fBm prior models in inverse or ill-conditioned problems. The following works should nevertheless be mentioned, although they consider prior properties of fBm, and in particular an empirical power-law spectrum decay, rather than considering prior fBm models themselves. In the context of restoration of remote sensing imagery [START_REF] Jalobeanu | Bayesian estimation of blur and noise in remote sensing imaging[END_REF] or of magnetic resonance imaging [START_REF] Tirosh | Polyharmonic smoothing splines for multi-dimensional signals with fractal-like spectra[END_REF], prior spectrum power-law decay have been used for the estimation of scalar fields from noisy or incomplete observations. Several recent works have also been led in this direction in the context of optic-flow estimation. In particular, prior power-law models for power spectra [START_REF] Derian | Wavelets to reconstruct turbulence multifractals from experimental image sequences[END_REF] or for structure functions [START_REF] Heas | Power laws and inverse motion modeling: application to turbulence measurements from satellite images[END_REF] have recently been proposed in the context of turbulent motion estimation from image sequences. This paper proposes a generic methodological framework for the regularization of ill-posed or ill-conditioned inverse problems using a bi-dimensional, and possibly divergence-free, fBm prior. The proposed model possesses a structure similar up to some extent to the regularizer proposed very recently in [START_REF] Tafti | On regularized reconstruction of vector fields[END_REF]. However, the modeling point of view adopted in this recent work is very different. Indeed, the authors do not intend to provide some prior model for fBm but rather propose regularization models satisfying invariance under coordinate transformations. Moreover, conversely to this related work, the present framework exploits an advantageous wavelet representation of fBm expressed in the spatial domain. Because of their advantageous whitening properties, divergence-free fractional Laplacian wavelet bases are chosen, yielding simple Gaussian white noise priors on wavelet coefficients. Furthermore, a low-complexity optimization algorithm, avoiding spectral computation and relying on mono-dimensional fast wavelet transforms, is proposed to maximize the posterior model.

The reminder of the paper is the following. After defining the properties of isotropic divergence-free fBm, we recall in the following section their spectral and wavelet representations. In section 3, we introduce a generic Bayesian framework for the regularization of inverse problems, based on the previous wavelet representation of fBm. Optimization algorithms for the maximization of the posterior model are then described in section 4. The modeling relevance and the algorithm efficiency are finally discussed in section 5, in the context of fractal image restoration and turbulent motion estimation.

2. Divergence-free isotropic fBm . For the clarity of the presentation, this work is restricted to the bi-dimensional case although it is in principle possible to extend the regularization to three dimensions. Moreover, in order to make the presentation as generic as possible, we focus on the more general case of divergence-free vectorial fields, keeping in mind that scalar fields or divergent vector fields are particular cases.

2.1. Divergence-free isotropic fBm definition and properties. The unique self-similar zero-mean Gaussian random process with stationary increments is fractional Brownian Motion (fBm) [START_REF] Flandrin | On the spectrum of fractional brownian motions[END_REF]. The definition of one-dimensional scalar fBm has been extended to the case of one-dimensional multivariate processes [START_REF] Amblard | Basic properties of the multivariate fractional Brownian motion[END_REF], to the case of isotropic -i.e., multi-dimensional-scalar fields [START_REF] Reed | Spectral representation of fractional brownian motion in n dimensions and its properties[END_REF], or to the case of isotropic multivariate signals [START_REF] Monin | Statistical Fluid Mechanics: Mechanics of Turbulence[END_REF][START_REF] Tafti | Fractional Brownian vector fields[END_REF]. We first recall the properties of the latter isotropic fBm model in the particular case of two-dimensional divergence-free vector fields.

Let u denote a realization1 of the bi-dimensional divergence-free isotropic fBm random variable u = (u 1 , u 2 ) T . Let us consider the Sobolev space of vectorial functions defined for ζ ∈ R by:

(H ζ (R 2 )) 2 = {f ∈ (L 2 (R 2 )) 2 : (1 -∆) ζ/2 f ∈ (L 2 (R 2 )) 2 }.
Realizations belong to the space of divergence-free functions:

u = (u 1 , u 2 ) T ∈ H ζ div (R 2 ) = {u ∈ (H ζ (R 2 )) 2 : div x u = 0},
where the divergence operator is defined in the sense of distribution. More precisely, for any test function ψ ∈ C 1 (R 2 ) with a fast decay at infinity and with vanishing moments, we define the divergence operator by:

div x u div x u, ψ L 2 (R 2 ) = -u/∇ x ψ , (2.1)
where ., . and ./. denote the inner product of two scalar (resp. two vectorial) functions in

L 2 (R 2 ) (resp. in (L 2 (R 2 )) 2
). Definition (2.1) is necessary since fBm realizations are not continuously differentiable [START_REF] Flandrin | On the spectrum of fractional brownian motions[END_REF]. Divergence-free isotropic fBm are vectorial random fields in H ζ div (R 2 ) with the following properties [START_REF] Tafti | Fractional Brownian vector fields[END_REF]: • zero-mean Gaussian variables, • non-stationary vector fields with homogeneous increments, • statistically invariant by rotation (isotropic),

• statistically self-similar with parameter H ∈ [0, 1].

Let us detail these properties. Statistical self-similarity and statistical rotation invariance are defined respectively by the following equality in probability law: for any orthogonal matrix M , i.e., M M T = I, ∀x ∈ R 2 :

u(σx) L =σ H u(x), ∀σ ∈ R+ (2.2) u(x) L =M u(M T x), (2.3)
where the real exponent H is known as the Hurst exponent.

The exponent has to be chosen as H ∈ [0, 1], since the Fourier integral defining the fBm process diverges for H > 1. Nevertheless, we will see that expanding fBm with wavelets with sufficiently vanishing moments enables the generalization to exponent

H ∈ R + /N [6].
Non-stationarity of the random vector field u is a direct consequence of its selfsimilarity [START_REF] Flandrin | On the spectrum of fractional brownian motions[END_REF]. Let E u [•] denote the expectation operator with respect to the probability distribution of u. The isotropic divergence-free fBm increments are homogeneous and characterized by the following auto-correlation function [START_REF] Tafti | Fractional Brownian vector fields[END_REF]:

E u [(u i (x + τ ) -u i (x -h + τ ))(u j (x) -u j (x -h)) * ] ∝ [G H (h + τ )] ij -2[G H (τ )] ij + [G H (h -τ )] ij , ∀H ∈ [0, 1], (2.4)
with ∝ representing the symbol equal up to a multiplicative constant and where [G H (h)] ij is the (i, j)-th component of the 2 by 2 matrix G H (h) defined by:

G H (h) = h 2H I - hh T h 2 .
(2.5)

The auto-correlation taken at τ =0 gives the second order moment of the increments:

E u [(u i (x) -u i (x -h))(u j (x) -u j (x -h)) * ] ∝ [G H (h)] ij , ∀H ∈ [0, 1], (2.6)
For H > 1, one can show that the corresponding Gaussian fields have stationary increments of order H + 1 where H denotes the largest integer lower or equal to H, see [START_REF] Tafti | Fractional Brownian vector fields[END_REF].

2.2. Spectral representation of fBm. Let dW (k) = (dW 1 (k), dW 2 (k)) T be the differential of a vectorial complex bi-dimensional Brownian motion satisfying the following properties:

dW i (-k) = dW * i (k), (2.7)
E W [dW i (k)] = 0, (2.8)
E W [dW i (k)dW * j (k)] = δ ij dk. ( 2
u(x) -u(0) = 1 (2π) 2 R 2 (e ik•x -1) k σG -H 2 (k)dW (k) (2.10)
where matrix G -H 2 is given by (2.5), where σ is a positive constant and dW (k) is the differential of a vectorial complex bi-dimensional Brownian motion satisfying (2.7)-(2.9). That is to say, assuming a vanishing initial condition u j (0) = 0, the j-th component writes:

u j (x) = 1 (2π) 2 R 2 (e ik•x -1)σ k -H-1 i (δ ij - k i k j k 2 )dW i (k). (2.11)
Proof. The construction of the fBm in the spectral domain is derived in an analogous manner as what has been done previously in the scalar isotropic case [START_REF] Reed | Spectral representation of fractional brownian motion in n dimensions and its properties[END_REF]. In order to be consistent with previous work [START_REF] Reed | Spectral representation of fractional brownian motion in n dimensions and its properties[END_REF][START_REF] Tafti | Fractional Brownian vector fields[END_REF], we first define the frequency representation of an isotropic divergence-free fractional Gaussian noise by the Wiener integral:

1 (2π) 2 R 2 e ik σG -H 2 (k)dW (k). (2.12)
where the so-called fractal frequency filter G -H 2 (k) is given by (2.5). Then, using the fact that fBm increment between two position-vector points x and y is given by the integration of the fractional noise along the straight-line path, connecting x and y (see [START_REF] Reed | Spectral representation of fractional brownian motion in n dimensions and its properties[END_REF]) we obtain: First, wavelets yield interesting fBm representations since non stationarity and self-similarity of fBm necessitate some appropriate space-scale analysis. Wavelets constitute whitening filters for these signals, which are correlated both in space and scale. The simplest choice is to use orthonormal wavelet bases. In this context, it can be shown that the variance of fBm wavelet coefficients decays through scales according to some power-law related to H [START_REF] Flandrin | Wavelet analysis and synthesis of fractional brownian motion[END_REF]. Nevertheless, with orthonormal basis, coefficients remain correlated in space and scale, and only approximate fBm whitening filters can be designed. A better choice is to use bi-orthogonal fractional wavelet bases. They represent ideal whitening filters [START_REF] Elliott | A fourier-wavelet monte carlo method for fractal random fields[END_REF][START_REF] Meyer | Wavelets, generalized white noise and fractional integration: the synthesis of fractional brownian motion[END_REF][START_REF] Tafti | Invariances, laplacian-like wavelet bases, and the whitening of fractal processes[END_REF][START_REF] Tafti | Fractional Brownian vector fields[END_REF]. Moreover, bi-orthogonal wavelet bases can also be designed to represent the space of divergence-free functions, which is a constraint of incompressible flows [START_REF] Lemarié-Rieusset | Analyses multirésolutions non orthogonales, commutation entre projecteurs et dérivation et ondelettes vecteurs à divergence nulle[END_REF].

u(x) -u(y) = 1 (2π) 2 R 2 (e ik•x -e ik•y ) k G -H 2 (k)dW (k). ( 2
Second, choosing wavelets to expand fBm removes the singularity at the zero frequency in (2.11) for any Hurst exponent H > 1 as long as we are considering wavelets with sufficient vanishing moments [START_REF] Tafti | Fractional Brownian vector fields[END_REF]. Convergence of the Fourier integral is therefore guaranteed, and the concept of fBm may be extended to "integrated" forms.

The third argument is related to the special context of optimization in Bayesian analysis: wavelet multiresolution analyses of unknowns proved experimentally to be relevant optimization strategy for non-convex inverse problems [START_REF] Dérian | Wavelet-based fluid motion estimation[END_REF]. Moreover, optimization can be led with a low-complexity since divergence-free wavelet decomposition [START_REF] Deriaz | Divergence-free and Curl-free wavelets in 2D and 3D, application to turbulence[END_REF] or the one-dimensional fractional wavelet decomposition [START_REF] Abry | The wavelet-based synthesis for fractional brownian motion proposed by f. sellan and y. meyer: Remarks and fast implementation[END_REF] have an advantageous fast implementation with recursive filter banks [START_REF] Mallat | A Wavelet Tour of Signal Processing: The Sparse Way[END_REF].

We postpone the latter point to section 4 and introduce in the following representations of isotropic bi-dimensional fBm with divergence-free and fractional Laplacian wavelets.

Construction with divergence-free wavelets.

To take into account implicitly the incompressibility constraint, the vectorial fBm u is represented in terms of divergence-free wavelet series [START_REF] Deriaz | Divergence-free and Curl-free wavelets in 2D and 3D, application to turbulence[END_REF]. The anisotropic 2 divergence-free wavelet spaces can be defined by:

H ζ div (R 2 ) = span{Ψ ,s ; , s ∈ Z 2 } (2.14)
with:

Ψ ,s = a s curl[χ ,s ] (2.15)
where curl

( ∂ ∂x2 , -∂ ∂x1 ) t , a s 1 √ 4 s 1 +4 s 2 and χ ,s = ψ 1,s1 ⊗ ψ 2,s2 represents an orthonormal scalar basis of L 2 (R 2 ), constructed using regular wavelets {ψ ,s = 2 s/2 ψ(2 s x-); , s ∈ Z}. Wavelets 3 Ψ ,s compose a primal vectorial basis of (L 2 (R 2 )) 2 .
Precisely, this primal divergence-free basis is constituted by wavelets

Ψ ,s = (Ψ 1 ,s , Ψ 2 ,s ) T , where Ψ 1 ,s = a s (ψ 1,s1 ⊗ ∂ ∂x2 ψ 2 ,s2 ) and Ψ 2 ,s = -a s ( ∂ ∂x1 ψ 1,s1 ⊗ ψ 2 ,s2
). Let us now consider measurements test functions constituted by the dual wavelets 4 Ψ ,s = 2 The anisotropy of the wavelet basis has not to be confused with the isotropy of the fBm model.

3 with the L 2 norm Ψ ,s 2 = R 2 i Ψ i ,s Ψ i ,s (x)dx = 1 4 with the L 2 norm Ψ ,s 2 = R 2 i Ψi ,s Ψi ,s (x)dx = 1 ( Ψ1 ,s , Ψ2 ,s
) T , which are biorthogonal to the primal wavelets Ψ ,s . Using such test functions, one can decompose uniquely any divergence-free vector field u ∈ H ζ div (R 2 ):

u = , s ∈ Z 2 u/ Ψ ,s Ψ ,s (x). (2.16)
Given this biorthogonal wavelet basis, the Wiener integral can be decomposed as an infinite sum of weighted independent standard Gaussian random variables, yielding the following result.

Proposition 2.2. Let u ∈ H ζ div (R 2
) be an isotropic fBm of parameter H, with ζ=-H-1 and let Ψ ,s be the Fourier transform of Ψ ,s . fBm u can be constructed by means of the divergence-free wavelet series:

u(x) = c ,s ,s R 2 (e ik•x -1) k ζ Ψ ,s (k)dk (2.17)
where c = σ/(2π) 2 is a positive constant and where the set of coefficients { ,s , ( , s) ∈ Z 2 } are distributed according to independent and identically distributed (i.i.d) zeromean Gaussians of variance equal to the L 2 norm of the dual wavelets, i.e., ,s ∼ N (0, Ψ ,s 2 ).

Proof. For any vectorial function

v = (v 1 , v 2 ) t ∈ span{ Ψ ,s ; , s ∈ Z 2 }
, since the wavelets Ψ ,s form a basis of this functional space, we have:

R 2 i vi (k)dW i (k) = R 2 ,s,i v/ Ψ ,s Ψi ,s (k)dW i (k) = ,s v/ Ψ ,s ,s (2.18)
with

,s = i R 2 Ψi * ,s (k)dW i (k), (2.19)
where Ψ ,s is the Fourier transform of Ψ ,s . This equality stands for any test function v. Therefore one obtains that the differential of Brownian motion can be decomposed as:

i dW i (k) = ,s,i ,s Ψi ,s (k)dk. (2.20)
Since the L 2 norm is conserved in the Fourier domain because of Plancherel's theorem, by definition of the Wiener integral (2.19) we get:

,s ∼ N (0, Ψ ,s 2
).

(2.21) Thus, using (2.16), (2.11) can be rewritten as:

u j (x) = 1 (2π) 2 ,s,i ,s R 2 (e ik•x -1)σ k -H-1 δ ji - k j k i k 2 Ψi ,s (k)dk. (2.22)
The term (δ ji -kj ki k 2 ) represents the symbol of Leray orthogonal projector, from (L 2 (R 2 )) 2 to the divergence-free function space H ζ div (R 2 ). Since divergence-free wavelets satisfy by construction i ,s (δ ji -

k j k i k 2 ) Ψi ,s = ,s Ψj ,s ,
this operator can be removed, yielding the divergence-free wavelet representation of isotropic fBm (2.17).

From proposition 2.2, we can derive the structure of the expected isotropic power spectrum, which is in agreement with [START_REF] Reed | Spectral representation of fractional brownian motion in n dimensions and its properties[END_REF].

Corollary 2.1. Let S(κ) represent the sphere around wave number κ. The mean expected two-dimensional power spectrum of u ∈ H ζ div (R 2 ), with ζ=-H-1 obeys to the following power-law:

E 2d (κ) 1 2π k S(κ= k ) E ,s [ c ,s ,s k -H-1 Ψ ,s (k) 2 ]dk = c 2 ,s E ,s [ ,s 2 ] k -2H-2 ∝ κ -2H-2 , (2.23)
resulting in the following expected isotropic power spectrum,

E(κ) 2πκE 2d (κ) ∝ κ -2H-1
When mentioned in the following, we may consider an isotropic power spectrum decay κ -2H-1 occurring on a bounded spectral interval Υ ⊂ R.

Construction with fractional Laplacian wavelets. The term k ζ Ψj

,s (k) in (2.17) represents in the Fourier domain the wavelet fractional Laplacian (-∆)

ζ 2 Ψ j
,s (x). Thus, (2.17) can be rewritten as a linear combination of fractional Laplacian wavelets, yielding a divergence-free fractional Laplacian representation:

u j (x) = c ,s ,s (-∆) ζ 2 Ψ j ,s (x) -Ψ j ,s (x) x=0 .
(2.24)

Indeed, one can define the fractional differentiation operator D ζ in the Fourier domain by [START_REF] Meyer | Wavelets, generalized white noise and fractional integration: the synthesis of fractional brownian motion[END_REF]:

Dζ = (ik) ζ (2.25)
and thus obtain the divergence fractional wavelets:

(-∆) ζ 2 Ψ j ,s (k) = (-(i k ) 2 ) ζ 2 Ψj ,s (k) = k ζ Ψj ,s (k), (2.26)
where (∆) ζ = (D 2 x + D 2 y ) ζ represents a fractional Laplacian differentation. Note that the fractional differentiation preserves the divergence-free property:

k • (-∆) ζ 2 Ψ ,s (k) = k • k ζ Ψ ,s (k) = k ζ k • Ψ ,s (k) = 0, (2.27)
with the notation (-∆)

ζ 2 Ψ ,s = a s curl[(-∆) ζ 2 χ ,s ]. Let us now introduce bi-dimensional vectorial functions of H ζ div (R 2 ): Φ ,si = (Φ 1 ,si , Φ 2 ,si ) T = 1 √ 2 curl[χ 0 ,si ] with χ 0 ,s1 = ψ 1,s1 ⊗ ϕ 2,0 and χ 0 ,s2 = ϕ 1,0 ⊗ ψ 2,s2
, where ϕ ,0 denote scaling functions at scale 2 -0 = 1 associated to the scalar orthonomal wavelets ψ ,s . According to the work of Meyer [START_REF] Meyer | Wavelets, generalized white noise and fractional integration: the synthesis of fractional brownian motion[END_REF], the low frequency contribution in the fBm representation of (2.24) can be isolated, yielding the following corollary to proposition 2.2:

Corollary 2.2. Let u ∈ H ζ div (R 2
) be an isotropic fBm of parameter H, with ζ=-H-1. fBm u can be constructed by means of the divergence-free fractional Laplacian wavelet series:

u j (x) = si>0, ,si [(-∆) ζ 2 Φ j ,si (x)] + c s1,s2>0, ,s [(-∆) ζ 2 Ψ j ,s (x)], (2.28)
where we have low-frequency function coefficents ,si , and where according to proposition 2.2 wavelet coefficients ,s are i.i.d zero-mean Gaussian variables.

3. Bayesian estimation.

3.1. Likelihood models. We will consider in this paper two inverse-problems: image restoration and divergence-free optic-flow estimation. Let us introduce likelihood models. Minus their logarithms correspond to so-called data-terms noted J l , which give a relation between observations and the variable of interest. Let y comprise either the observed image {y 0 } to be enhanced or the image pair {y 0 , y 1 } used for optic-flow estimation. In their simplest form, the restoration data-term writes:

J l (u, y) = 1 2 R 2 (u(x) -y 0 (x)) 2 dx, u ∈ L 2 (R 2 ), (3.1)
and, for optic-flow estimation, the Displaced Frame Difference (DFD) data-term is defined as:

J l (u, y) = 1 2 R 2 (y 1 (x + u(x)δt) -y 0 (x)) 2 dx, u ∈ (L 2 (R 2 )) 2 (3.2)
In order to simplify notations we will consider in the following the case δt = 1. The maximum likelihood estimate can be defined as the functional minimum, i.e., arg min u J l (u, y), subject to the divergence-free constraints div x u(x) = 0, ∀x ∈ R 2 in the optic-flow case. However, it is well-known that both problems are ill-posed and need to be regularized.

3.2. Prior models . Using (2.21), simple regularization models follows from decomposition (2.28): an isotropic divergence-free fBm prior can simply be formalized by independent, identically distributed, zero-mean and normalized Gaussian distributions for fractional wavelet coefficients:

u j (x) = si>0, ,si [(-∆) ζ 2 Φ j ,si (x)] + c s1,s2>0, ,s [(-∆) ζ 2 Ψ j ,s (x)], Λ ∼ N (0, Λ), (3.3)
where denotes a vector composed by the set of coefficients { ,s , ,s } and Λ is a diagonal matrix with unitary entries for coefficients ,s whose scale are related to the spectral range Υ, and with zero entries otherwhere. The zero entries on the diagonal of Λ implies that the multivariate Gaussian prior has an infinite variance in these directions. Note that the probabilistic model (3.3) implies the power-law decay of power spectrum E(κ) ∝ κ -2H-1 , ∀κ ∈ Υ.

3.3. Posterior estimation. The computation of unknown u is now reduced to the estimation of its fractional divergence-free wavelet coefficients, knowing that a subset of the coefficients ,s are realization of independent one-dimensional reduced and zero-mean normal distributions.

MAP estimate. An interesting estimate may be obtain solving the minimization problem:

* = arg min J l ( , y), where Λ ∼ N (0, Λ) (

with

J l ( , y) = 1 2 R 2 (g • z (x) -y 0 (x)) 2 dx,
where in the case of an optic-flow problem we define (z (x) where J p , so-called regularization term, is defined as minus the logarithm of the prior, i.e.,

-x) ∈ H ζ div (R 2 ): z (x) = x + si>0, ,si [(-∆) ζ 2 Φ j ,si (x)] + c s1,s2>0, ,s [(-∆) ζ 2 Ψ ,s (x)] with g(x) = y 1 (x), or in the case of an image restoration problem z ∈ H ζ (R 2 ): z : (x) → si>0, ,si [(-∆) ζ 2 χ j ,si (x)] + c s1,s2>0, ,s [(-∆) ζ 2 χ ,s (x)
J p ( ) = 1 2 T Λ . (3.7)
Note that by a change of variable, (3.5) is solved equivalently by minimization with respect to independent and identically distributed Gaussian variables c of constant variance 1 βc 2 .

3.3.3. On equivalence with other regularizers. The following proposition clarifies the equivalence of the regularizer (3.7) with the model proposed very recently in [START_REF] Tafti | On regularized reconstruction of vector fields[END_REF]. Furthermore, it shows that a particular case of (3.7) is the vorticity derivative penalization, which constitutes a popular approach in the context of fluid flow estimation [START_REF] Corpetti | Dense estimation of fluid flows[END_REF][START_REF] Yuan | Discrete orthogonal decomposition and variational fluid flow estimation[END_REF]. Let the operator curl † denote the formal adjoint of the scalar curl operator. It is given for any vectorial function

v = (v 1 , v 2 ) t ∈ (L 2 (R 2 )) 2 , by the vectorial product: curl † v = (∂ x1 v 2 -∂ x2 v 1 )
. The scalar quantity curl † v is called in the following the vorticity of the vector field v.

Proposition 3.1. Let u ∈ H ζ div (R 2 )
. For an identity matrix Λ, regularization with (3.7) is equivalent to penalization of the L 2 norm of the signal fractional Laplacian derivative:

J p ( ) = 1 2 (-∆) -ζ 2 u 2 (3.8)
In the particular case ζ = -2, we have:

J p ( ) = 1 2 ∇curl † u 2 = ∆(u) 2 .
(3.9)

A consequence of this proposition is that, in two-dimension, an appropriate regularization for an isotropic divergence-free (resp. scalar) fBm of parameter H = 2 (resp. H = 1) is the penalization of the L 2 norm of the vorticity gradient (resp. of the Laplacian).

Proof. We have:

(-∆) -ζ 2 u 2 R 2 j |(-∆) -ζ 2 u j (x)| 2 dx = R 2 j | , s ∈ Z 2 u/(-∆) ζ 2 Ψ ,s (-∆) -ζ 2 Ψ j ,s (x)| 2 dx = , s ∈ Z 2 u/(-∆) ζ 2 Ψ ,s 2
T , which proves (3.8) and the equivalence with the regularizers proposed in [START_REF] Tafti | On regularized reconstruction of vector fields[END_REF]. Now, let us highlight the fact that penalization of the L 2 norm of the Laplacian or of the vorticity gradient are equivalent in the two-dimensional divergence-free case. Indeed, let χ be a stream function such as:

u = curl χ. (3.10)
From the definition of vorticity in 2D, one can prove that:

-∆χ = curl † (u), (3.11) which implies that:

∇curl † (u) 2 = -∇∆χ 2 = ∆(curl χ) 2 = ∆(u) 2 .
(3.12)

4. Optimization. Optimization strategies analogous to the algorithm proposed by [START_REF] Dérian | Wavelet-based fluid motion estimation[END_REF] are used to solve the minimization problem (3.6). Optimization relies on a wavelet multi-resolution analysis based at each level on a LBFGS procedure, i.e., on a quasi-Newton gradient descent algorithm subject to the strong Wolf conditions [START_REF] Nocedal | Numerical Optimization[END_REF]. The regularization term is quadratic. Thus, in theory, for Gaussian likelihood models such as (3.1), the algorithm converges to the global minimum, while for non-Gaussian models such as (3.2), it is only insured to converge to the nearest local minima. Interesting wavelet-based multiresolution strategies dealing with these non-linearities have been recently proposed in [START_REF] Dérian | Wavelet-based fluid motion estimation[END_REF]. The gradient descent optimization requires at each step the evaluation of the functional to be minimized and of its gradient.

We hereafter propose two alternative optimization strategies, which either imply spectral or spatial computation. The spectral optimization method, relying on an exact fractional Laplacian, is efficient for low order fractional derivatives. However in practice, for high order derivatives (typically for ζ > -1), fractional wavelet bases become poorly conditioned, resulting in inefficient optimization algorithms. Nevertheless, we propose for such situations (typically for ζ < -1) an efficient spatial optimization method, relying on a fractional Laplacian approximation. For the clarity of the presentation, we will restrict in the following to optic-flow functionals although it is easy to derive analogous expressions in the case of image restoration problems.

Spectral optimization method.

4.1.1. Functional analytical gradient. The following proposition proves that the functional gradient can simply be obtained by Fast Wavelet Transform (FWT).

Proposition 4.1. The gradient with respect to vector of the functional minimized in (3.6) is the vector composed of wavelet coefficients, given by FWT of the two-dimensional scalar function

βa s curl † [c (-∆) ζ 2 (g • z -y 0 ) (∇ x g) • z ](x)
using the orthonormal wavelet basis {χ ,s , ( , s) ∈ Z 2 }, supplemented by the additional vector Λ .

Proof. The gradient of the data-term J l ( , y) is given (in the case of the divergencefree optic-flow problem) by inner products with the fractional divergence-free wavelets:

∂ ,s J l ( , y) = c [(g • z -y 0 ) (∇ x g) • z ](x)/a s curl[(-∆) ζ 2 χ ,s ](x) ∂ ,s i J l ( , y) = c [(g • z -y 0 ) (∇ x g) • z ](x)/a s curl[(-∆) ζ 2 χ ,si ](x) . (4.1)
Using Parseval's theorem we get:

∂ ,s J l ( , y) = c [ (g • z -y 0 ) (∇ x g) • z ](k)/a s i(k 2 , -k 1 ) t |k| ζ χ ,s (k)] = c |k| ζ a s i(-k 2 , k 1 )[ (g • z -y 0 ) (∇ x g) • z ](k), χ ,s ](k) = c a s curl † [(-∆) ζ 2 (g • z -y 0 ) (∇ x g) • z ](x), χ ,s ](x) . (4.2)
The second inner product of (4.1) can be computed analogously. Moreover, it is straightforward to show that the gradient of the regularization term is simply:

∂ J p ( ) = Λ . (4.3)
4.1.2. Spectral domain computation. Based on (4.2), we derive a spectral method for the computation of the gradient of functional minimized in (3.6). It is based on Fast Fourier Transform (FFT) and FWT:

Algorithm 1. (functional gradient) i) FFT of the components of [(g • z -y 0 ) (∇ x g) • z ](x),
ii) compute in Fourier domain the fractional differentiation and apply the curl † operator on

[(g • z -y 0 ) (∇ x g) • z ](k) iii) inverse FFT to get curl † [(-∆) ζ 2 (g • z -y 0 ) (∇ x g) • z ](x)
iv) decompose the result by FWT using the orthogonal wavelets a s χ ,s and the scaling functions χ 0 ,si / √ 2 to obtain the data-term gradient (4.2), v) functional gradient obtained by adding vector (4.3)

In order to evaluate the functional or its gradient, one needs to reconstruct the unknown u appearing in g • z or (∇ x g) • z from the fractional divergence-free wavelet coefficients . This can be done by spectral computation. Indeed, (2.28) can be rewritten as:

u = 1 (2π) 2 (R 2 ) 2 (dk, dk) T e ik•x i(k 2 , -k 1 ) T |k| ζ (4.4) R 2 dx( si>0, ,si √ 2 χ 0 ,s +c s1,s2>0, ,s a s χ ,s )e -ik•x ,
yielding the following reconstruction algorithm:

Algorithm 2. (reconstruction of fBm) i) inverse FWT of using the orthogonal wavelets a s χ ,s and scaling functions 1 √ 2 χ 0 ,s , ii) FFT of this scalar function iii) fractional differentiation and curl operator in Fourier domain, iv) inverse FFT of the components to obtain u Algorithms 1 and 2 yield the ingredients necessary to approach the MAP estimate * with a gradient descent method.

Spatial optimization method.

In order to avoid spectral computation of fractional differentiation, i.e., preserve the low-order complexity, and by the same way, prevent numerical instabilities for large negative decay ζ, we choose to rotate the posterior in a standard divergence-free wavelet basis.

Rotation of the posterior. Let d denotes the set of divergence-free wavelet coefficients {d

,s = u/ Ψ ,s = a -1 s χ, χ ,s ; ( , s) ∈ Z 2 }.
With these notations, the reconstruction formula (2.16) becomes for the vector field u and its potential χ:

u = , s ∈ Z 2 d ,s Ψ ,s (x); χ = , s ∈ Z 2 d ,s a s χ ,s (x).
(4.5)

The MAP estimation problem (3.6) is therefore rewritten as a minimization problem with respect to wavelet coefficients d:

d * = arg min d {βJ l (d, y) + J p (d)} (4.6) with J p (d) = 1 2 (d) T Λ (d), (4.7)
Therefore, fractional divergence-free wavelet coefficients need to be expressed in terms of standard divergence-free wavelet coefficients d. ,s , ,s .Twodimensional divergence-free fractional Laplacian wavelet coefficients are approached by a linear combination of divergence-free wavelet coefficients, implying low-dimensional matrices products :

[ ] = im i=0 -α i R (-α-i) T d a -1 s R (i) , i m ∈ N * (4.8)
where [ ] and [d γs ] represent respectively the matrix of two-dimensional fractional coefficients and two-dimensional standard coefficients d ,s multiplied by factor γ s .

Proof. Since the potential χ and fBm u have the same wavelet decomposition coefficients up to the factor a s , we can write that:

,s = u/(-∆) ζ 2 Ψ ,s = a -1 s ,s d ,s χ ,s , (-∆) ζ 2 χ ,s (4.9)
where (-∆)

ζ 2 χ ,s denotes one dual component of the fractional scalar wavelet basis.

Although operator (-∆)

ζ 2 is not separable, bi-dimensional connection coefficients appearing in (4.9) can be approached by generalized Newton's binomial expansion5 :

χ ,s , (-∆) ζ 2 χ ,s = ∞ i=0 -ζ 2 i r (-ζ 2 -i) 1 ,s 1 , 1,s1 r (i) 2 ,s 2 , 2,s2 (4.10)
with the mono-dimensional connection coefficients:

r (-ζ 2 -i) 1 ,s 1 , 1 ,s1 = ψ 1 ,s 1 , -∂ 2 ∂x 2 -ζ 2 -i ψ 1 ,s1 , (4.11)
r (i) 2 ,s 2 , 2 ,s2 = ψ 2 ,s 2 , -∂ 2 ∂x 2 i ψ 2,s2 . (4.12)
See appendix for details on wavelet connection coefficients computation. The infinite sum in (4.10) needs in practice to be truncated. Let us choose connection coefficients vanishing for large translation, i.e., for | 1 -1 | → ∞ (see appendix for more details), corresponding to the following truncation order:

i m = max{i ∈ N|(-ζ -2i + 1) ∈ R + }.
In practice, as illustrated in figures of appendix A, for sufficiently negative exponent ζ, binomial coefficients with i > i max weighting the expansion (4.10) are negligible, i.e., the proposed truncation yields good approximations. Now, using the separable truncated expansion (4.10), (4.9) can be rewritten as low-dimensional matrices products.

We have:

,s ≈ a -1 s im i=0 -ζ 2 i 1 ,s 1 r (-ζ 2 -i) 1 ,s 1 , 1,s1 2 ,s 2 d 1 , 2 ,s 1 ,s 2 r (i) 2 ,s 2 , 2 ,s2 (4.13) [ ] ≈ im i=0 -ζ 2 i R (-ζ 2 -i) T d a -1 s R (i) .
Therefore, the MAP estimation problem (4.6)-(4.7) becomes:

d * = arg min d {βJ l (d, y) + J p (d)} (4.14) with J p (d) = 1 2 I Λ { im i=0 -ζ 2 i R (-ζ 2 -i) T d a -1 s R (i) } 2 F ,
where A 2 F = trace(A T A) is the Frobenius norm of matrix A and with

I Λ {[ ]} 2 F = t Λ .
4.2.2. Functional analytical gradient. We derive hereafter a trackable expression of the gradient of functional minimized in (4.14). Proposition 4.3. Let us define (in the case of optic-flow problems) the functions

z d (x) -x ∈ H ζ div (R 2 ) : x → ,s d ,s Ψ ,s (x
) and g(x) = y 1 (x). The gradient with respect to vector d of the functional minimized in (4.14) is the vector composed of wavelet coefficients, given by FWT of the two-dimensional vectorial function

βc g • z d -y 0 (∇ x g) • z d
using the dual divergence-free wavelet basis { Ψ ,s , ( , s) ∈ Z 2 }, supplemented by the additional vector whose ( , s)-th component writes

  im,jm i,j=0 -ζ 2 i -ζ 2 j R (-ζ 2 -i) I Λ {R (-ζ 2 -j) T d a -2 s R (j) }R (i) T   ,s
.

Proof. Analogously to results in [START_REF] Dérian | Wavelet-based fluid motion estimation[END_REF], it is straightforward to see that the gradient of the data-term J l (d, y) is given by inner-products with divergence-free wavelets:

∂ d ,s J l (d, y) = c [ g • z d -y 0 (∇ x g) • z d ](x)/Ψ ,s (x) . (4.15)
As a matter of fact, inner products in (4.15), i.e., the data-term gradient, are simply given by a FWT of c g • z d -y 0 (∇ x g) • z d using the dual basis of divergence-free wavelets. Note that an analogous expression can be obtained for image restoration defining

z d : (x) ∈ L(R 2 ) → c ,s d ,s χ ,s (x) ∈ L(R 2
). Concerning the regularization term, using the matrix equality:

∂ [d] 1 2 Trace(R (i) T [d] T R (-ζ 2 -i) R (-ζ 2 -j) T [d] R (j) T ) = R (-ζ 2 -i) R (-ζ 2 -j) T [d]R (j) R (i) T , (4.16)
it follows from the linearity of the trace operator that the gradient of the regularization term writes:

∂ d ,s J p (d) = im,jm i,j=0 -ζ 2 i -ζ 2 j R (-ζ 2 -i) I Λ {R (-ζ 2 -j) T [d a -2 s ]R (j) }R (i) T ,s , (4.17)
where (A) ,s denotes the ( , s)-th element of matrix A.

Spatial domain computation.

Using proposition 4.3, a trackable computation of the gradient of functional minimized in (4.14) is provided by means of simple FWT and low-dimensional matrix products. Moreover, reconstruction (2.16) by inverse FWT is straightforward for any values of coefficients d. We have therefore provided the two ingredients, which are necessary to approach the MAP estimate d * with a gradient descent algorithm.

Experiments.

5.1. Divergence-free isotropic fBm generator. Realizations of isotropic fBm were generated using the Fourier-based reconstruction formula (4.4), yielding scalar or divergence-free vectorial fields. More precisely, in agreement with the fBm model (3.3), the wavelets coefficients { ,s } were sampled according to centered and reduced Gaussian white noise. Coefficents { ,s }, corresponding to the low frequency contribution, were neglected and arbitrarily fixed to zero. Nevertheless, in principle, they should have been drawn according to a fractional ARIMA random walk , see [START_REF] Meyer | Wavelets, generalized white noise and fractional integration: the synthesis of fractional brownian motion[END_REF]. Isotropic fBm were then synthesized from this Gaussian white noise by application of algorithm 2.

The three following fBm realizations were synthesized in order to form an evaluation benchmark for the regularization model: fBm-1. An isotropic fBm scalar field of 128 × 128 pixels, synthesized using this procedure, is displayed in figure 5.1. Pixel grey levels have been normalized. The wavelet generator was constructed from orthogonal Daubechies with 4 vanishing moments (Daubechie-4) with periodic boundary conditions. The Hurst exponent was fixed to the limit case of H = 0, implying in the fBm model (3.3) a fractional Laplacian wavelets (-∆) ). This powerlaw spectrum is visible in figure 5.5 (resp. figure 5.8). These synthesized fields are representative of self-similar incompressible turbulent motion in the 2D and 3D case, as described by the Kolmogorov theory [START_REF]The local structure of turbulence in inompressible viscous fluid for very large reynolds number[END_REF]. • Block Matching 3D (BM3D). This algorithm is very popular in computer vision. It has been recently proposed in [START_REF] Dabov | Image denoising by sparse 3-d transformdomain collaborative filtering[END_REF]. It is based on nonlocal patchwise image modeling and achieves outstanding performance in the restoration of natural scenes. The true noise variance was provided as an input to the algorithm. • Penalization of L 2 norm of Laplacian. We used the Matlab implementation of this well-known deblurring algorithm. The optimal regularization parameter was chosen (using a brute-force approach). Power spectrum of fBm-1 with H = 0 (pink dashed curve), spectrum of the noisy observation (red crosses) and power law decay in E(κ) ∝ κ -1 (yellow dashed line). Estimated power spectra using BM3D algorithm (turquoise squares), Laplacian L 2 norm (green crosses) or fractional Laplacian (blue stars).

• Penalization of L 2 norm of fractional Laplacian (3.8) with ζ = -1. We solved (3.6) using the proposed Fourier-based optimization (see algorithm 1 and algorithm 2), since it is efficient for such a low order fractional derivative (conversely to the spatial optimization procedure, see appendix A for details). A wavelet generator identical to the one used for fBm-1 generation was employed, i.e., Daubechie-4 wavelets with periodic boundary conditions. Not surprisingly, the true prior achieved the best performances in terms of PSNR as shown in the following table.

PSNR

Obs. BM3D Laplacian Frac. Laplacian fBm. 1 23.3957 26.1707 26.4211 26.4846 Moreover, as illustrated in figure 5.1, the fractional Laplacian prior succeeded to restitute the true regularity of fBm-1, while other methods seemed to fail. As shown in figure 5.2, the diagnostic is confirmed in the Fourier domain: BM3D algorithm and the Laplacian regularization over-estimate the isotropic power spectrum decay while under-estimating energy at small scales. On the contrary, the fractional Laplacian model reconstruct the true power spectrum with a good accuracy in the entire scale range.

In summary, the previous experiments highlight the fact that, efficient image restoration methods for computer vision scenes, are not necessarily adapted to fractal image restoration. An appropriate modeling is required for these irregular stochastic processes. 2) were generated by deformation of an image y 0 with fBm-2 or fBm-3, i.e., by the following interpolation y 1 (x) = y 0 (x -u(x)). We used cubic B-splines for this effect. Pixel grey levels of each of the two image couple (y 0 , y 1 ) were then normalized. The input images y 0 and y 1 are displayed in figure 5.3 for fBm-2 (resp. figure 5.6 for fBm-3), together with their associated ground truth motion and vorticity.

The vectorial divergence-free fields fBm-2 and fBm-3 were estimated solving the minimization problem (4.6), combining the DFD data-term (3.2) with four different regularizers. Therefore, for all of these experiments, estimation was carried out in terms of coefficients of the reconstructed motion divergence-free wavelet series. The wavelet generator was constructed from divergence-free biorthogonal Coiflet-10 with periodic boundary conditions. A non negligible benefit of using wavelets is their capability to provide an efficient multiresolution analysis framework, which enables the recovery of large displacements, as shown in [START_REF] Dérian | Wavelet-based fluid motion estimation[END_REF]. The four different regularization models are listed hereafter:

• High-order polynomial approximation. An implicit regularization by polynomial approximation is a well-known approach in computer vision. We used here local approximations by polynomials of degree 10 based on coiflet vanishing moments, as proposed in [START_REF] Dérian | Wavelet-based fluid motion estimation[END_REF]. • Penalization of L 2 norm of velocity components gradients. The most common approach in optic-flow estimation is to penalize the velocity components gradients L 2 norm as first proposed in [START_REF] Horn | Determining optical flow[END_REF]. We used the wavelet-based implementation proposed in [START_REF] Kadri-Harouna | Divergence-free wavelets and high order regularization[END_REF] with an optimally tuned regularization coefficient (using a brute-force approach).

• Penalization of L 2 norm of vorticity gradient. In fluid motion estimation a popular approach is to penalize the L 2 norm of the vorticity gradient, i.e., to use model (3.9) [START_REF] Corpetti | Dense estimation of fluid flows[END_REF][START_REF] Yuan | Discrete orthogonal decomposition and variational fluid flow estimation[END_REF]. However, as remarked previously, this prior is the true regularization model for fBm only in the particular case when H = 0. This is not the case neither for fBm-2 nor for fBm-3. We used the waveletbased implementation proposed in [START_REF] Kadri-Harouna | Divergence-free wavelets and high order regularization[END_REF] with an optimally tuned regularization coefficient (using a brute-force approach). • Penalization of L 2 norm of fractional Laplacian (3.8) with ζ = -3 for fBm-2 (resp. ζ = -4/3 for fBm-3). We used the spatial optimization method (see proposition 4.3) since it is efficient for those high order Laplacian fractional Power spectrum of divergence-free isotropic fBm-3 (turquoise dashed curve) and power law decay in E(κ) ∝ κ -5/3 (yellow dashed line). Estimated power spectra using polynomial approximation (red crosses), velocity gradient penalization (green crosses), vorticity gradient penalization (blue stars), fractional Laplacian penalization (pink squares).

derivative (see appendix A for justifications).

Estimated fBm fields were evaluated based on the Root Mean Squared end-point Error (RMSE) and the Mean Barron Angular Error (MBAE), see [START_REF] Baker | A database and evaluation methodology for optical flow[END_REF] Using appropriate fBm priors, that is to say fractional Laplacian regularizers, yields clearly the best results in terms of RMSE or MBAE criteria. The enhancement in reconstruction accuracy brought by those priors is more striking than in the previous image restoration experiments. The choice of the regularization model is likely to have here a greater influence, since optic-flow estimation is not a closed problem conversely to image restoration. Consequently, if some physical knowledge about the signal is available, e.g., if we know that the flow can be represented by Kolmogorov's turbulence model [START_REF]The local structure of turbulence in inompressible viscous fluid for very large reynolds number[END_REF], it is crucial to design accordingly some appropriate fBm prior. We now compare and focus on vorticity fields computed from estimated motion fields. Indeed, the quality of the reconstruction of this quantity is an important matter in turbulence. Figure 5.4 (resp. figure 5.7) shows that the vorticity regularity of fBm-2 (resp. of fBm-3) is better preserved using the proposed regularization. In particular, these figures show the limitation of standard approaches: penalizing the gradient of velocity or the gradient of vorticity yields coherent vorticity blobs with smooth edges, which are not representative of fBm. As shown in figure 5.5 and figure 5.8, inspection of power spectra computed with inappropriate priors reveals an over-estimation (resp. an under-estimation) of the energy of fBm-2 (resp. of fBm-3) at intermediate scales (resp. at small scales), i.e., a non respect of the spectrum power-law decay. Note that for the estimation of fBm-2, the slight deviation induced by the proposed regularization from the true spectrum at small scales has a negligible influence in terms of RMSE and MBAE since the spectrum decay is very important in this case.

Conclusion.

In this paper we a prior model for regularization of isotropic divergence-free fBm of parameter H. Realizations of these stochastic processes belong to the Sobolev function space (H ζ (R 2 )) 2 , with ζ=-H-1 and satisfy weak divergence-free property. We show that these probabilistic models can be constructed by means of divergence-free fractional Laplacian wavelet series and prove that a subset of coefficients of the series are distributed according to Gaussian white noise. We then propose to use these prior models for regularization of ill-conditioned problems in a Bayesian framework. We point out some particular cases where the proposed model corresponds to well-known or more recent regularizers. We design two optimization strategies in order to access the MAP estimate by gradient descent algorithms. The first method relies on the computation by FFT and FWT of fractional Laplacian wavelets and is efficient for small values of H. It exploits the very simple expression of fBm prior with fractional Laplacian wavelet series. We complement this Fourierbased method by a low complexity optimization algorithm which is efficient for large values of H. It only implies FWT computation and low-dimensional matrix products. It exploits the expression of these fractional priors rotated in standard divergence-free wavelet bases by the mean of fractional Laplacian wavelet connection coefficients. We finally show weaknesses of state-of-the-art regularizers and the usefulness of these fBm prior models in the context of fractal image restoration and estimation of turbulent motion fields from image couples. Using the framework proposed by Beylkin, we obtain with this kernel the recursion [START_REF] Beylkin | On the representation of operator in bases of compactly supported wavelets[END_REF][5]:

r ( α 2 ) ,0 = 2 α L-1 k=0 k-L-1 j=k h[k]h[k -j]r ( α 2 )
j+2 ,0 , (A.4) subject to the asymptotics of r

( α 2 )
,0 for large :

r ( α 2 ) ,0 =        1 c α 1+α + O 1 1+α+2n for > 0, 1 c α (-) 1+α + O 1 (-) 1+α+2n
for < 0.

(A.5)

By a direct inversion procedure, we solve linear system (A.4), which is of infinite dimension, but where coefficients have been approached by the asymptotic laws (A.5) for sufficiently large translations. Note that, for large negative translations, the asymptotic of connection coefficients differs from the asymptotic of standard fractional differentiation coefficients of [START_REF] Beylkin | On the representation of operator in bases of compactly supported wavelets[END_REF].

1 Z exp -1 2 T

 12 ] with the identity function g(x) = x. From a Bayesian perspective, (3.4) corresponds to the following Maximum A Posteriori (MAP) problem: * =arg max p(y| )p( ) (3.5) where p(y| ) = 1 Z y| exp -βJ l ( ,y) , p( ) = Λ with parameter β ∈ R + * , and where Z y| (resp. Z ) denote the partition function (independent of ) of the likelihood p(y| ) (resp. of the centered, of identical variance and uncorrelated Gaussian prior p( )) distribution. 3.3.2. Regularized functional. The solution of (3.5) can simply be obtained by minimizing minus the logarithm of the posterior distribution: p( |y) ∝ p(y| )p( ), yielding the MAP estimate: * = arg min {βJ l ( , y) + J p ( )} (3.6)

Proposition 4 . 2 .

 42 Let R (α) = [r (α) ,s , ,s ], α ∈ R/Z -denote the matrix of monodimensional fractional Laplacian wavelet connection coefficients r (α) ,s , ,s = ψ ,s , composed at line indices ( , s ) and column indices ( , s) by the elements r (α)

ζ 2

 2 with ζ = -1. Consequently, it induces in (2.23) a power spectrum behavior E(κ) ∝ κ -1 in a scale range going from the pixel size to the image dimension, i.e., Υ = [1/128, 1]. This power-law spectrum is visible in figure 5.2. fBm-2 and fBm-3. Two divergence-free vectorial fBm fields of 128×128 for fBm-2 (resp. 256 × 256 for fBm-3) pixels were also synthesized using this procedure. They are displayed in figure 5.3 and figure 5.6. The wavelet generator was constructed from divergence-free biorthogonal Coiflet with 10 vanishing moments (Coiflet-10) and periodic boundary conditions. The Hurst exponent was fixed to H = 2 (resp. H = 1/3), implying in the fBm model a fractional Laplacian wavelets (-∆) ζ 2 with ζ = -3 (resp. ζ = -4/3). Consequently, it induces in (2.23) a power spectrum behavior E(κ) ∝ κ -5 (resp. E(κ) ∝ κ -5/3 ) in a scale range going from the pixel size to the image dimension, i.e., Υ = [1/32, 1/16] (resp. Υ = [1/256, 1]

Fig. 5 . 1 .

 51 Fig. 5.1. Above: fBm-1 realization u(x) (left) and observed noisy fBm-1 y(x) (right). Below from left to right: estimates obtain by penalizing Laplacian L 2 norm, using BM3D algorithm and penalizing fractional Laplacian L 2 norm.

5. 2 .

 2 Fractal image restoration. An additive zero mean i.i.d Gaussian noise was added to the synthesized scalar field fBm-1. It resulted in a Peak Signal to Noise Ratio (PSNR) of 23.3957 dB. The scalar fBm was restored using the data-term (3.1) combined with three different regularizers/algorithms :

  Fig.5.2. Power spectrum of fBm-1 with H = 0 (pink dashed curve), spectrum of the noisy observation (red crosses) and power law decay in E(κ) ∝ κ -1 (yellow dashed line). Estimated power spectra using BM3D algorithm (turquoise squares), Laplacian L 2 norm (green crosses) or fractional Laplacian (blue stars).

Fig. 5 . 3 .

 53 Fig. 5.3. From left to right. Images y 0 (x) and y 1 (x), synthesized divergence-free isotropic fBm-2, vorticity of fBm-2.

Fig. 5 . 4 .

 54 Fig. 5.4. fBm-2 vorticity MAP estimation (from left to right) polynomial approximation, velocity gradient penalization, vorticity gradient penalization and fractional Laplacian penalization.

5. 3 .

 3 Turbulent divergence-free optic-flow estimation. Images y 1 in the data-term (3.

Fig. 5 . 5 .

 55 Fig. 5.5.Power spectrum of divergence-free isotropic fBm-2 (turquoise dashed curve) and power law decay in E(κ) ∝ κ -5 (yellow dashed line). Estimated power spectra using polynomial approximation (red crosses), velocity gradient penalization (green crosses), vorticity gradient penalization (blue stars), fractional Laplacian penalization (pink squares).

Fig. 5 . 6 .

 56 Fig. 5.6. From left to right. Images y 0 (x) and y 1 (x), synthesized divergence-free isotropic fBm-3, vorticity of fBm-3.

Fig. 5 . 7 .

 57 Fig. 5.7. fBm-3 vorticity MAP estimation using (from left to right) polynomial approximation, velocity gradient penalization, vorticity gradient penalization and fractional Laplacian penalization.

Fig. A. 1 .

 1 Fig. A.1. Above: binomial coefficients in (A.1) for ζ=-4 (red dots indicate non-neglected coefficients, i.e., i ≤ im=3). Below: connection coefficients r (-ζ 2 -i) ,0 (left) and r (-i) ,0 (right) in function of , for i ≤ im (from top to bottom i=1, 2 and 3).

Fig. A. 2 .

 2 Fig. A.2. Above: binomial coefficients in (A.1) for ζ=-7/3 (red dots indicate non-neglected coefficients, i.e., i ≤ im=2). Below: connection coefficients r (-ζ 2 -i) ,0 (left) and r (-i) ,0 (right) in function of , for i ≤ im (from top to bottom i=1 and 2).

Fig. A. 3 .

 3 Fig. A.3. Above: binomial coefficients in (A.1) for ζ=-1 (red dots indicate non-neglected coefficients, i.e., i ≤ im=1). Below: connection coefficients r (-ζ 2 -i) ,0 (left) and r (-i) ,0 (right) in function of , for i ≤ im (i = 1).

2 )

 2 coefficients. Indeed, (4.13) can be rewritten as: elements are connection coefficients between standard ϕ(x) and fractional Laplacian scaling functions at the finest scale:Let us note that the truncated series (A.1) constitutes a good approximation for sufficiently large negative values of exponents ζ. Indeed as shown in figure A.1, A.2 and A.3, relation (A.1) is exact for ζ = -4 or almost exact for ζ = -7/3, while it constitutes a poor approximation for ζ = -1. For α/2 ∈ N, fractional Laplacian operator becomes a standard differentiation, up to factor (-1) α/2 . In order to easily compute (A.2) in a more general case, the fractional Laplacian operator is rewritten as a convolution operator. Indeed, if (α + 1)/2 ∈ R/N, fractional Laplacian can also be defined by Riesz potential6 [15]: 1+α)/2) . In the previous expression, the convolution kernel writes k(x) = 1 c α |x| α+1 .

  . Results are reported in the following table.

	RMSE/MBAE	Polynomial	Gradient	Vorticity gradient Frac. Laplacian.
	fBm.2	0.706/19.543 0.563/14.759	0.530/14.314	0.454/11.859
	fBm.3	0.213/ 7.159	0.195/7.266	0.212/7.234	0.156/5.305

Realizations of random variable are noted with calligraphic symbol. For example, a realization of the random variable u is noted u.

Generalized Newton's binomial expansion stands for any non-negative integer: -ζ/2 / ∈ Z -. However, this is never the case since in(2.22) we have set -ζ/2 = H + 2, with H > 0.

This definition can be extend to the case (α + 1)/2 ∈ N using some appropriate kernel, see e.g.,[START_REF] Reichel | Characterization of balls by riesz-potentials[END_REF]

Acknowledgments. The authors would like to thank Frédéric Lavancier for his helpful comments, which have helped to improve the paper.

Appendix A. Fractional Laplacian wavelet connection coefficients. Note that by simple bi-dimensional wavelet transform, one can derive the matrices of wavelet connection coefficients from the matrix of scaling function connection