Henri-Alex Esbelin
email: alex.esbelin@univ-bpclermont.fr

Rémy Malgouyres
email: remy.malgouyres@udamail.fr

Sparse Convolution-based Digital Derivatives, Fast Estimation For Noisy Signals and Approximation Results

Keywords: Digital Analysis, Convolution, Digital Derivative, Derivative Estimation, Curvature Estimation

We provide a general notion of a Digital Derivative in 1-dimensional grids, which has real or integer-only versions. From any such masks, a family of masks called skipping masks are defined. We prove general results of multigrid convergence for skipping masks. We propose a few examples of digital derivative masks, including the now well-known binomial mask. The corresponding skipping masks automatically have multigrid convergence properties. We study the cases of parametric curves tangents and curvature. We propose a novel interpretation of digital convolutions as computing points on a smooth curve, the regularity of which depends on the mask. We establish, in the case of binomial and B-spline masks, a close relationship between the derivatives of the smooth curve, and the digital derivatives provided by the mask.

Introduction

Aim of the work

The goal of this work is to put in place a foundation stone for a new paradigm for numerical calculations. The motivation for this is that the objects of current methods of numerical analysis, which deal mathematically primarily with real numbers, cannot effectively be computed, or even stored in memory with state of the art computers, even taking account the possibilities of more powerful machines implementing Turing computation in a foreseeable future. This problem has been tackled in several ways.

The most obvious, most widely used way is floating point computations. The numbers are encoded by approximations which lie in a (non translation invariant) set of, so called, representable numbers. When any operator is applied, be it an arithmetical operation or a function (say an analytical function) of a number, the result, which is most often not representable, is approximated, in a way which is not always well specified in norms and may be implementationdependent. During the execution of a complex numerical method, accumulated rounding-off errors can create convergence and unstability problems for the method, even if the mathematical method is proved convergent, and stable with respect to initial values. But most of all, it is generally very difficult to prove that the result is correct with such methods.

Interval Analysis (see [START_REF] Alefeld | Interval analysis: theory and applications[END_REF] for instance) consists in performing numerical computation by applying operators to interval, rather than to (encoding of) numbers. Each operation on an interval consists of operations, or possibly bounds calculation which are representable and calculable, for the boundaries of the interval. In this way, the technique allows to obtain bounds for solutions of equations or numerical problems. In many cases, however, the estimations are too loose to be really useful.

Exact Real arithmetic (see an introduction in [START_REF] Vuillemin | Exact real computer arithmetic with continued fractions[END_REF]) rests on the notion of a computable number, which consists in restricting constructions of real numbers by completion of the set Q of rational numbers (by Cauchy Sequences, or by decreasing sequences of intervals), narrowing the construction to calculable (through some Turing Complete language) sequences. We can define an arithmetic on the sets of numbers which, under some restrictions related to computability, implement exactly numerical computation on calculable real numbers.

Symbolic Calculus consists in formal manipulations algebraic expressions, by rewriting formally the expressions according to the rules of the algebraic structures involved. This approach to exact real computation is generally considered a natural complement to exact real arithmetic, which can fill in gaps leaved by non-computable operations. This approach, however, is quite remote from numerical calculus as we purpose to do.

In our approach, we would define (hopefully) any numerical problem in a discrete form and in a continuous form. The discrete form of the problems would involve only integers (with some constant denominator or scale factor). Each of the considered operator would possibly change the rational scale factor (which can be stored), but otherwise be computable using only integers. Then, a multigrid convergence result would allow us to show that the integeronly computable solution of the discrete problem converges to the real valued solution of the continuous (and original) version of the problem. This avoids the accumulation of rounding-off errors typical of floating point computations.

To make such an approach succesfull, we would need to redefine so called digital, or discrete versions of all the objects involved in real analysis. Having made a short wish list (quadratic optimisation, gradient method, PDE integration, and so on), we found that the notion of a derivative is at the root of all these notions. In this paper, we propose a general definition of a digital derivative, and prove some multigrid convergence results.

In the framework of image and signal processing, as well as shape analysis, a common problem is to estimate derivatives of functions, or tangents and curvatures of curves and surfaces, when only some (possibly noisy) sampling of the function or curve is available from aquisition. This problem has been investigated through finite difference methods, scale-space ([KSW05], [START_REF] Lindeberg | Scale-space for discrete signals[END_REF]), and digital geometry ( [START_REF] Coeurjolly | Multigrid convergent principal curvature estimators in digital geometry[END_REF], [START_REF] Coeurjolly | Integral based curvature estimators in digital geometry[END_REF], [START_REF] Kerautret | Curvature estimation along noisy digital contours by approximate global optimization[END_REF], [START_REF] Lachaud | Analysis and comparative evaluation of discrete tangent estimators[END_REF], [START_REF] Lachaud | Fast, accurate and convergent tangent estimation on digital contours[END_REF], [START_REF] Levallois | Parameter-free and multigrid convergent digital curvature estimators[END_REF], [START_REF] Mazo | About multigrid convergence of some length estimators[END_REF], [START_REF] Provot | Estimation of the derivatives of a digital function with a convergent bounded error[END_REF], [START_REF] De | Comparison and improvement of tangent estimators on digital curves[END_REF]).

In our previous work ( [START_REF] Esbelin | Convergence of binomial-based derivative estimation for c 2 noisy discretized curves[END_REF], [START_REF] Esbelin | Taylor optimal kernel for derivative etimation[END_REF], [START_REF] Malgouyres | Binomial convolutions and derivatives estimation from noisy discretizations[END_REF], [START_REF] Fourey | Normals estimation for digital surfaces based on convolutions[END_REF], [START_REF] Gonzalez | Convergence of level-wise convolution differential estimators[END_REF]) an approach to derivative estimation from digitized data was introduced. As in scale-space approaches, this approach is based on simple computations of convolutions. However, unlike scale-space methods, this approach is oriented towards integer-only models and algorithms, and is based on a discrete approach to analysis on Z. Implementation of the convolution-based approach is straightforward. As far as the speed of convergence is concerned, the method for tangents of [START_REF] Malgouyres | Binomial convolutions and derivatives estimation from noisy discretizations[END_REF] was proved in [START_REF] Esbelin | Convergence of binomial-based derivative estimation for c 2 noisy discretized curves[END_REF] uniform worst-case O(h 2/3 ) for C 3 functions, where h is the size of the pixel. Moreover, our estimator allowed some (uniformly bounded or stochastic) noise. Furthermore, the method allowed to have a convergent estimation of higher order derivatives, and in particular a uniform O(h 4/9 ) estimation for the curvature of a generic curve.

The weakness of these results was the high complexity of computation. In this paper, we prove that the same convergence rates are allowed with a very low complexity by using sparse masks. This idea was previously used in connection with Taylor optimal masks ([EM14]). It turns out that it is more efficient with sparse Binomial masks.

Outline

The paper is organized as follows. In Section 2, we present the notion of a digitization of a realvalued function, in which a discrete sequence of sample values approximates a real function. We also present a few noise models, including basic quantization, uniform noise or bias (depending on the scale), and stochastic noise.

In Section 3, we define digital derivatives. This notion is based on digital convolutions with masks which satisfy some identities involving sums (so-called moments). We also define some important characteristics, such as the convergence order of the mask. We also show some good properties, notably related to composition of digital derivatives by convolving the corresponding masks.

In Section 4, we begin to investigate errors for the purpose of proving multigrid convergence results. We distinguish the sampling error, which is due to the fact that we know only the values of the functions at a (locally) finite set of samples, from the values error, which is due to inaccurate input values (due to quantization, noise or bias). We also give upper bounds for both types of errors, for all the considered types of noise models.

In Section 5, we study some smoothing masks, i.e. some 0-derivative masks. These masks are important because, basically, any 1-derivative mask can be seen as a smoothing mask composed with a usual simple finite difference mask. The intuitive idea is that if we use the simple finite difference mask on noisy data, the importance of noise will be very detrimental to the results. By smoothing, which is nothing but a weighed average of values of the function, some of the noise cancels out (at least for a stochastic noise with expected value 0). We present three families of masks which we consider important, among the many masks on which we worked until now ([EMC11], [START_REF] Esbelin | Taylor optimal kernel for derivative etimation[END_REF], [START_REF] Malgouyres | Binomial convolutions and derivatives estimation from noisy discretizations[END_REF], and [START_REF] Gonzalez | Convergence of level-wise convolution differential estimators[END_REF] in which we began to seek low computational complexity for convolutions). These three families of masks are the Binomial, Taylor-Optimal, and B-spline masks.

In Section 6, we establish multigrid convergence results of the digital derivatives toward the corresponding continuous derivatives of the original function, providing in addition a worst case rate of convergence. That section contains the main convergence results of this paper, namely uniform convergence results in the uniform noise or bias model, and convergence results of the standard deviation of our estimation in the stochastic model. Both convergence results, which hold for derivatives of any order ω, have a rate O(h α-ωα/k 0 ), where 0 < α ≤ 1 depends on the noise and k 0 ≥ 2 depends on the digital derivative operator. Furthermore, these convergence results apply to a general family of masks called skipping masks, which are basically obtained from any ω-derivative mask, by using a small mask with few values (for better computational complexity), but using a sparse set of sample values at points of a regular sub-grid.

In Section 7, we define three families (also called the Binomial, Taylor-Optimal, and B-spline masks) of skipping 1-derivative masks, with as corollary (from the results of Section 6) some multigrid convergence results obtained directly simply by checking the definition of a digital 1-derivative.

In Section 8.2, we apply these operators and results to some additional particular cases, such as second derivative, or the estimation of tangents and curvatures of parametrized curves.

In Section 9, we present a novel interpretation of digital convolutions as the computation of some points on a smooth curve, the differentiability of which depends on the kernel. We study in details the relationship between binomial kernels for derivatives of any order and the corresponding derivatives of Bézier curves. We do the same for so-called B-spline kernels and B-spline curves. Since the B-spline case is somewhat simpler, we present it first. The B-spline kernels, which we introduce in this paper for the first time, seem a bit ad-hoc. We believe, however, that the B-spline case, just as the Bézier case, will be of great theoretical interest in the future.

The last section is devoted to experiments, the main purpose of which are manly to illustrate the difference between masks, but also that only a few (typically a small constant such as 10 to 20 points) values are needed to obtain a good estimation. An appendix also proposes different criteria for choosing a mask, depending on the considered application. Functions for which domain and range are sub-ring of R, without any assumption on the nature of these sets, are called real functions. We call digital function a function from Z to Z.

Contents

First we fix a model for the relationship between a continuous function and its digitization. Let f : R -→ R be a real function and let Γ : Z -→ Z be a digital function. Let h be the digitization step (i.e. the size of a pixel). We introduce a (possibly noisy) digitization of f . Definition 2.1 The function Γ is a digitization of f with error ϵ h with digitization step h on the interval I if for any integer i such that ih ∈ I we have:

hΓ(i) = f (ih) + ϵ h (i) (1)
In order to apply this theory to signal analysis, we could distinguish two digitization steps for the domain and for the range, using the following definition:

h 2 Γ(i) = f (ih 1 ) + ϵ h 1 ,h 2 (i).
The results concerning the estimation errors bounds may be easily generalized to this framework. However generalized convergence results may hardly be efficient because of the many parameters involved.

We consider the following particular models for the errors ϵ h on the values:

• Exact Values In this model, the values are known exactly:

ϵ h = 0
Note that, although this model has been the most widely used in approximation theory (e.g. [START_REF] Michael | On the convergence of derivatives of bernstein approximation[END_REF], [START_REF] Michael | On the deviation of a parametric cubic spline interpolant from its data polygon[END_REF]), this value error model is not very realistic from an Information Sciences point of view.

• Uniform Noise (or Uniform Bias) on Values: In this model, the error ϵ h on the values is uniformly bounded by some constant which depends on the sampling step h. In our model, however, this bound can be asymptotically greater that h. Namely we assume here that

0 ≤ |ϵ h (i)| ≤ Kh α
where 0 < α ≤ 1 and K is a positive constant. Note that this error can also have some bias, in the sense that the average noise value (or expected value) could be non-zero. This model has been considered by authors working on digital estimators based on thick digital lines.

• Quantization of Values: In this model, the errors ϵ h on the values is uniformly bounded by 1 2 h. This is a particular case of uniform noise with α = 1, and corresponds to the case when some basic quantization has been obtained by rounding-off the exactly known values of the function, for example for digital storage. This case is equivalent to

Γ(i) = [ f (ih) h ] .
A variant is when quantization has been obtained by an integer part (floor case):

0 ≤ ϵ h (i) < h , which is equivalent to Γ(i) = ⌊ f (ih) h ⌋.
This model is the most usual for digital estimators based on digital lines.

• Stochastic Noise on Values: In this model, the errors ϵ h (i) on the different values are independent random variables with expected value 0 and standard deviation σ(h), converging to 0 along with h. For convergence results, we shall suppose σ(h) ≤ Kh α . This model, and generally some estimators using measure theory, have been used recently ( [START_REF] Cuel | Robust geometry estimation using the generalized voronoi covariance measure[END_REF]), but, to the best of our knowledge, without any multigrid convergence results.

Digital Derivatives for Functions

We introduce now a notion of digital derivatives.

Definition 3.1 (Digital Derivative) A digital ω-derivative mask is a sequence u = (u(i)) i∈Z ∈ R Z of finite support such that ∑ i∈Z i k u(i) = 0 for 0 ≤ k < ω and ∑ i∈Z i ω u(i) = (-1) ω ω! Its associated digital ω-derivative operator is the function ∆ u with domain R Z and co-domain R Z defined by ∆ u (v)(n) = (u ⋆ v)(n) = ∑ i∈Z u(i)v(n -i).
Its convergence order is the greatest integer ρ such that for all integers k such that 1 + ω ≤ k ≤ ρ, we have

∑ i∈Z i k u(i) = 0.
We shall sometimes say a derivative mask as a shorthand for a digital derivative mask.

Remark 3.1 Note that if there is no integer k such that 1 + ω ≤ k and

∑ i∈Z i k u(i) = 0, then the convergence order ρ is equal to ω.
Remark 3.2 In the case when the values u(i) are rational numbers and the sequence v has integer values, then the values ∆ u (v)(n) have rational values with denominator less than or equal to the GCD of the denominators of the u(i).

Remark 3.3 Note that this definition of a derivative operator easily extends to some families of sequences with infinite support, under some convergence hypothesis (e.g. rapidly decreasing masks like a sampling of the Gaussian). So, this definition can be extended to include scalespace notions [START_REF] Lindeberg | Scale-space for discrete signals[END_REF]).

For the sake of clarity, we shall denote by ∆ ω u the associated digital derivation operator of u when u is a digital ω-derivative mask. Intuitively, and as appears in the proofs below, the order ρ of an ω-derivative operator is the number of coefficients of the Taylor development of the function

1 h ω-1 ∆ ω u (f )-f (ω)
which vanish, minus one. As shown below, the convergence order ρ will determine the speed of convergence of ω) , which is to distinguish from the derivative order ω. Notice that the complexity of the computation is mainly dependent on the non zero coefficients of the mask, hence on the cardinality of its support (as long as it is computable). In [START_REF] Esbelin | Convergence of binomial-based derivative estimation for c 2 noisy discretized curves[END_REF], this cardinal was an increasing function of h. One novelty of this work is considering masks with supports of constant cardinality, which is in practice so small as 10 or 12. Hence the computation needs no more than that number of products of integers.

1 h ω-1 ∆ ω u (f ) to f (
We shall use the following generalization of the classical finite difference operator with order 1 and skipping step l: Definition 3.2 The following operator ∆ + l is a digital 1-derivative operator:

∆ + l (v)(n) = 1 l (v(n + l) -v(n))
the mask of which is defined by ∆ + l (-l) = 1 l , and ∆ + l (0) = -1 l and ∆ + l (i) = 0 for others values of i. Similarly, we denote by ∆ - l the following derivative operator:

∆ - l (v)(n) = 1 l (v(n) -v(n -l))
As we could expect, the composition of two digital derivatives is also a digital derivative:

Proposition 3.1 Let u = (u(i)) i∈Z be an ω-derivative mask and v = (v(i)) i∈Z be an ω ′ - derivative mask. Then u ⋆ v is a ω + ω ′ -derivative mask. Proof. Let k ≤ ω + ω ′ . Then ∑ n∈Z n k (u ⋆ v)(n) = ∑ n∈Z (i + (n -i)) k ∑ i∈Z u(i)v(n -i) = ∑ n∈Z ∑ i∈Z ∑ k p=0 ( k p ) i p (n -i) k-p u(i)v(n -i) = ∑ k p=0 ( k p ) (∑ n∈Z ∑ i∈Z i p u(i)(n -i) k-p v(n -i) ) = ∑ k p=0 ( k p ) ∑ j∈Z ∑ i∈Z i p u(i)j k-p v(j) = ∑ k p=0 ( k p ) ( ∑ j∈Z j k-p v(j) ) (∑ i∈Z i p u(i) )
This is zero except if k = ω + ω ′ and in this case all the terms are zero except if p = ω and in this case the sum is

( ω+ω ′ ω ′ ) (-1) ω ω!(-1) ω ′ ω ′ ! = (-1) ω+ω ′ (ω + ω ′ )! 2
The following remark can be easily checked by induction on the degree of a polynom:

Remark 3.4 Our definition of a digital ω-derivative corresponds partially to the usual derivative when applied to polynomial sequences: let u be a digital ω-derivative mask with convergence order ρ and let v(i) be defined by i k for a given non negative integer k; then for k ≤ ρ, we have

( ∆ ω u ((i k ) i∈Z ) ) (n) = k(k -1)...(k -ω + 1)i k-ω .
4 Basic Error Decomposition and Upper Bounds

Errors Related to Sampling and to Input Values

In order to show that the digital ω-derivative of a digitization Γ of a real function f provides an estimate for the continuous derivative f (ω) of f , we would like to evaluate, at each sample point, the difference between the digital derivative

1 h ω-1 (∆ ω u ⋆ Γ)(n)
of the digitized signal and the value of the usual ω-derivative f (ω) (nh) of f . This difference may easily be decomposed from Equation (1) and Definition 3.1 as the sum

1 h ω-1 (∆ ω u ⋆ Γ)(n) -f (ω) (nh) = ES ω (f, h, Γ, u, n) + EV ω (f, h, Γ, u, n) (2)
where

ES ω (f, h, Γ, u, n) = ( 1 h ω ∑ i∈Z u(i)f ((n -i)h) ) -f (ω) (nh) (3)
is called the sampling error, and

EV ω (f, h, Γ, u, n) = 1 h ω ∑ i∈Z u(i)ϵ h (n -i) (4)
is called the values error. As their names imply, the sampling error is due to the fact that we only know about the values of f at some grid points, and the values error is due to the fact that we do not know the exact values of f at sample points. The sampling error is a real values sequence. Under the uniform bias hypothesis, the values error is also a real valued sequence, but under the stochastic hypothesis, the values error is a sequence of random variable.

Upper Bound for the Sampling Error

In the following lemma, we show that the sampling error can be bounded independently from the error on input values, using the mask values, the norm of the k 0 -derivative of f , and a power of the digitization step (namely h k 0 -ω ), where k 0 is (at least in the case when f is sufficiently regular) the convergence order of the mask u. The consequence is some convergence results in the case when exact values of the function at sample points are known.

Lemma 4.1 Let f be C k on R, with k ≥ 1 + ω. Let u be a digital ω-derivative mask with convergence order ρ. Let k 0 = inf {k, 1 + ρ}. Let Γ be a digitization of f with step h on R. Suppose that f (k 0 ) is bounded on R. Then for all n, |ES ω (f, h, Γ, u, n)| ≤ h k 0 -ω ∥f (k 0 ) ∥ ∞ k 0 ! ( ∑ i∈Z |i k 0 u(i)| ) Proof.
From Taylor formula, we first give an upper bound for the sampling error:

ES ω (f, h, Γ, u, n) = 1 h ω ∑ i∈Z u(i) ( i=k 0 -1 ∑ j=0 (-ih) j j! f (j) (nh) + (-ih) k 0 k 0 ! f (k 0 ) (ζ(i)) ) -f (ω) (nh)
where ζ(i) lies between nh and (n -i)h. From the definition of the order, we get:

ES ω (f, h, Γ, u, n) = h k 0 -ω (-1) k 0 k 0 ! ( ∑ i∈Z i k 0 u(i)f (k 0 ) (ζ(i))
)

|ES ω (f, h, Γ, u, n)| ≤ h k 0 -ω ∥f (k 0 ) ∥ ∞ k 0 ! ( ∑ i∈Z |i k 0 u(i)| ) 2 Example 4.1 From Definition 3.1, if f is C 2 and a 1-derivative mask u has convergence order 1, then the sampling error is O(h). This is the case for the usual finite difference (∆(f )) (n) = f (n + 1) -f (n).
Example 4.2 Suppose that f is C 3 and that a 1-derivative mask u has convergence order 2, then the sampling error is O(h 2 ). This is the case for the usual symmetric finite difference (∆(f )) (n) = 1 2 (f (n + 1) -f (n -1)). For a 2-derivative mask v with convergence order 2, the sampling error for the second derivative is O(h).

Upper Bound for the Values Error

The following lemma gives an upper bound for the error related to uniform noise or uniform bias on the values at sample points (see Section 2). Lemma 4.2 Let f be a real function. Let u be a digital ω-derivative mask with convergence order ρ. Let Γ be a digitization of f with step h on R and uniform bias

ϵ h such that ∥ϵ h ∥ ∞ ≤ Kh α with 0 < α ≤ 1. Then for all n, |EV ω (f, h, Γ, u, n)| ≤ K h ω-α ( ∑ i∈Z |u(i)| ) Proof.
We straightforwardly give an upper bound for the values error:

|EV ω (f, h, Γ, u, n)| ≤ ∥ϵ h ∥ ∞ h ω ∑ i∈Z |u(i)| ≤ K h ω-α ( ∑ i∈Z |u(i)| ) . 2
The following lemma gives an upper bound for the error related to statistic noise with expected values 0 on the values at sample points.

Lemma 4.3 Let f be derivable on R. Let u be a digital ω-derivative mask. Let Γ be a digitization of f with step h on R and stochastic noise on values ϵ h with expected value 0 and standard deviation σ(h). Then for all n, the random variable

1 h ω-1 (∆ ω u ⋆ Γ)(n) -f (ω) (nh) has expected value ES ω (f, h, Γ, u, n) and standard deviation σ(h) h ω √∑ i∈Z (u(i)) 2 .
In other words, and roughly speaking, the global error is in this case statistically close to the sampling error.

Proof.

1 h ω-1 (∆ ω u ⋆ Γ)(n) -f (ω) (nh) is equal to the sum of the constant value ES ω (f, h, Γ, u, n) and the random variable defined by EV ω (f, h, Γ, u, n) = 1 h ω ∑ i∈Z u i ϵ h (n -i).
From linearity, its expected value is zero. Since the random variables ϵ h (n -i) are supposed to be independent and the support of u is supposed to be finite, the square of the standard deviation of EV ω (f, h, Γ, u, n) is equal to the sum of the squares of the standard deviations of

u i h ω ϵ h (n -i) which are |u i | h ω σ(h).
2 Note that for a fixed mask, the values error (or its standard deviation) generally does not converge to zero when h converges to 0. We shall propose below a way to make it tend to zero by adapting the mask to the digitization step (see Theorem 6.1 and Theorem 6.2 below).

Smoothing masks

A smoothing mask is a 0-derivative mask. In this case, the values (∆ 0 u ⋆Γ)(n) may be considered as partially de-noised values of f (nh). More precisely, under the uniform hypothesis, according to our upper bounds from Lemmas 4.2 and 4.1, the difference between h.(∆ 0 u ⋆ Γ)(n) and f (nh) cannot be considered not better (in order of magnitude) than the difference between h.Γ(n) and f (nh). Moreover, Lemma 4.3 with ω = 0 shows that, in the case of a stochastic noise on values, h.(∆ 0 u ⋆ Γ)(n) -f (nh) is a random variable, the expected value of which converges to 0 along with h, and its standard deviation is less or equal than the one of h.Γ(n) -f (nh).

Various choices for the smoothing masks may be considered. Here we focus on three families of masks.

Binomial Smoothing Masks

Definition 5.1 (Binomial Smoothing Masks) Let m be a positive integer and i be an integer. The binomial smoothing mask with diameter (2m -1) and cardinality 2m is the finite sequence of dyadic rational numbers A unified definition may be B (0) for r ≥ 0, are discrete smoothing masks with convergence order 0 if r is odd, and with convergence order 1 if r is even.

B (0) 2m-1 defined by B (0) 2m-1 (i) = 1 2 2m-1 ( 2m-1 m-i ) if -(m -1) ≤ i ≤ m,
(i) = 1 2 2m ( 2m m-i ) if -m ≤ i ≤ m,
r (i) = 1 2 r ( r ⌈ r 2 ⌉-i ) if -⌊ r 2 ⌋ ≤ i ≤ ⌊ r 2 ⌋,
Proof. We shall use the two following classical identities:

i=n ∑ i=0 i ( n i ) = n2 n-1 and i=n ∑ i=0 i 2 ( n i ) = n(n + 1)2 n-2 . ∑ i∈Z B (0) 2m-1 (i) = ∑ i∈Z B (0) 2m (i) = 1 are well known. ∑ i∈Z iB (0) 2m (i) = 0 come from parity. ∑ i∈Z iB (0) 2m-1 (i) = 1 2 2m-1 m ∑ i=-m+1 i ( 2m -1 m -i ) = 1 2 2m-1 m ∑ j=0 (m -j) ( 2m -1 j ) ∑ i∈Z iB (0) 2m-1 (i) = m - 1 2 (2m -1) = 1 2 ∑ i∈Z i 2 B (0) 2m (i) = 1 2 2m m ∑ i=-m i 2 ( 2m m -i ) = 1 2 2m m ∑ j=0 (m -j) 2 ( 2m j ) ∑ i∈Z i 2 B (0) 2m (i) = m 2 -2m 2 + 2m(2m + 1) 4 = m 2 2

Taylor-Optimal Smoothing Masks

The following so called Taylor-optimal smoothing mask is chosen having the best order as possible for a given diameter:

Definition 5.2 (Taylor-Optimal Smoothing Masks) Let m be a positive integer and i be an integer. The Taylor-optimal smoothing mask with diameter 2m and cardinality 2m is the finite sequence of rational numbers

O (0) 2m defined by O (0) 2m (i) = (-1) i+1 ( m |i| ) ( m+|i| m ) = (-1) i+1 ( 2m m+|i| ) ( 2m m ) if -m ≤ i ≤ m ,
2m is a digital smoothing mask with convergence order 2m -1.

Proof. Let k ≥ 0. Let f be a polynomial and define the classical finite difference operator

δ + by δ + (f ) (x) = f (x + 1) -f (x).
Iterating this operator leads to

δ 2m + (f ) (x) = 2m ∑ i=0 ( 2m i ) (-1) 2m-i f (x + i)
Applying this equality to the polynomials defined by

(x -n) k for 0 ≤ k ≤ 2m -1, we get 2m ∑ i=0 ( 2m i ) (-1) i (i + x -m) k = 0
Let now x = 0, then for all 0 ≤ k ≤ 2m -1 we have

0 = 2m ∑ i=0 ( 2m i ) (-1) i (i -m) k = (-1) m j=m ∑ j=-m ( 2m m + j ) (-1) j j k
Hence for k = 1, we have

∑ j=m j=-m,j̸ =0 (-1) j ( 2m m+j ) = - ( 2m m )
and for 2 ≤ k ≤ 2m -1, we have ∑ j=m j=-m,j̸ =0 (-1) j j k ( 2m m+j ) = 0. The same computation for the polynomial defined by (x -m) 2m leads to ∑ j=m j=-m,j̸ =0 (-1) j j 2m ( 2m m+j ) = (-1) m ̸ = 0. 2

B-splines Smoothing Masks

Now, we define the B-spline Smoothing Kernels.

Definition 5.3 Let s ∈ Z. The Cox-De-Boor functions N s,r (with knots t s = s), for r ∈ N, are defined inductively as follows.

N s,0 (t) = { 1 if s ≤ t < s + 1 0 otherwise.
For r ≥ 1, we use the induction formula (called the Cox-de Boor formula):

N s,r (t) = t -s r N s,r-1 (t) + s + r + 1 -t r N s+1,r-1 (t) Lemma 5.1 We have N i,r (t) = 0 for t ̸ ∈]i, i + r + 1[.
Besides, for all t ∈ R and r ∈ N, we have the following properties: Partition of Unity Property:

∑ i∈Z N i,r (t) = 1 Shift Property N s,r (t + 1) = N s-1,r (t) Symmetry Property N s,r ( s + r+1 2 + t ) = N 0,r ( s + r+1 2 -t )
Proof. These properties are well known and follow immediately by induction from the Cox-de Boor formula. 2 Definition 5.4 (B-spline Smoothing Masks) Let us consider a positive integer r ≥ 0.

We introduce a convolution mask (S (0) r (i)) i∈Z with support included in {-⌊ r+1 2 ⌋, . . . , ⌊ r+1 2 ⌋}, based on the Cox-de Boor functions with degree r on R, with (bi-infinite) nodal vector τ = (i) i∈Z . The mask values are defined by S (0)

r (i) = N 0,r ( r+1 2 + i) = N -i,r ( r+1 2 ) if -⌊ r+1 2 ⌋ ≤ i ≤ ⌊ r+1 2 ⌋ and S (0)
r (i) = 0 otherwise. This mask is called the uniform B-spline smoothing mask with degree r.

Since N 0,r ( r+1 2 ) = 0, the cardinality of the mask S (0)

r is 2⌊ r+1 2 ⌋.
The mask is always symmetric.

Proposition 5.3 The mask S (0)

r is a 0-derivative mask (i.e. a smoothing mask) with convergence order 1.

Proof. The fact that S (0) r is a 0-derivative mask follows immediately from the partition of Unity property in Lemma 5.1 and Definition 3.1. The convergence order (same definition) equal to 1 follows from:

∑ i∈Z i.S (0) r (i) = ∑ i∈Z i.N 0,r ( r + 1 2 + i )
The last sum is equal to zero from the symmetry property in Lemma 5.1. 2

The reader can find in Appendix A.3 some practical estimation for the constant involves in the asymptotic order sampling error for the three families of smoothing masks presented above.

Skipping Masks: Cheap Multigrid Convergence

The idea is adapting the mask to the step of digitization, in order to get 1

h ω-α (∑ i∈Z |u(i)|
) converging to zero along with h. For limiting the complexity of computation, we set the number of non zero coefficients of the mask fixed. This idea was independently used in [START_REF] Mazo | About multigrid convergence of some length estimators[END_REF] for length estimation. Other attempts have been considered in [START_REF] Gonzalez | Convergence of level-wise convolution differential estimators[END_REF].

Definition 6.1 (Skipping Masks) Let u be an ω-derivative mask with convergence order ρ, and l ∈ N * be a skipping step. The corresponding ω-derivative l-skipping mask u l is defined by u l (i) = 1 l ω u( i l ) if l divides i and 0 in any other cases. l is the length of the skips. As ∑ i∈Z i k u l (i) = l k-ω ∑ j∈Z j k u(j), the mask u l is an ω-derivative mask with same convergence order than u. But as

∑ i∈Z |u l (i)| = 1 l ω ∑ j∈Z |u(j)|
, it allows a convenient choice of l depending on h in order to get a zero converging values error.

Uniform Multigrid Convergence with Uniform Noise or Bias

Theorem 6.1 Let u be an ω-derivative mask with order ρ and u l the corresponding ω-derivative l-skipping mask with skips of length l. Suppose that f : R -→ R is a C k function and

f (k) is bounded, α ∈]0, 1], K ∈ R * + and h ∈ R * + . Suppose Γ : Z -→ Z is such that |hΓ(i) -f (hi)| ≤ Kh α for all i. Let k 0 = M in{k, 1 + ρ}. Then for k ≥ 1 + ω and l(h) = ⌊ h -1+α/k 0 ⌋ , we have ( 1 h ω-1 ∆ u l(h) ⋆ Γ ) (n) -f (ω) (nh) ∈ O(h α-ωα/k 0 ),
and, for k = ω and l such that lim h→0 h.l(h) = 0 and lim h→0 h α ((h.l(h)) ω = 0, we have

( 1 h ω-1 ∆ u l(h) ⋆ Γ ) (n) -f (ω) (nh) ∈ o(1).
Proof. First, we give an upper bound for the values error. From Lemma 4.2 and definitions, we have

EV (f, h, Γ, u l(h) , n) ≤ K l(h) ω h ω-α ∑ j∈Z |u(j)|. If l(h) = ⌊ h -1+α/k 0 ⌋ , it is easy to verify that 1 l(h) ω h ω-α ≤ h α-ωα/k 0 1 -h ω-ωα/k 0 which is O(h α-ωα/k 0 ).
We now turn to the sampling error:

• Let k ≥ 1 + ω. From Lemma 4.1 and definitions, the sampling error is less than

(l(h)h) k 0 -ω ∥f (k 0 ) ∥∞ k 0 ! ∑ j∈Z |j k 0 u(j)|. If l(h) = ⌊ h -1+α/k 0 ⌋ , it is easy to verify that |ES(f, h, Γ, u l(h) , n)| ≤ h α-ωα/k 0 ∥f (k 0 ) ∥∞ k 0 ! ∑ j∈Z |j k 0 u(j)|. • Let f be C ω on R. From Taylor formula, ES(f, h, Γ, u l(h) , n) turn to be equal to ∑ j∈Z j ω u(j)ϵ((n -jp(h))h), where ϵ is a function such that lim x→0 ϵ(x) = 0. 2
The multiplicative constants in the first part of Theorem 6.1 depends on various parameters. Some depend on informations about the data: a bound h max of the digitization step, or ∥f (k 0 ) ∥ ∞ or K depending on the model of the noise. Some others depend on the mask:

∑ i∈Z |i k u i | is named the k-moment
of the mask. In the case when k ≥ 1 + ω, the error is bounded by

( K 1 -h α-ωα/k 0 max ∑ i∈Z |u j | + ∥f (k 0 ) ∥ ∞ k 0 ! ∑ i∈Z |j k 0 u j | ) h α-ωα/k 0
For practical purposes, it could be useful to choose the mask in order to minimize this constant. In Appendix A.2, we provide formulas for some of our masks, and table of values for others.

Stochastic Multigrid Convergence with Stochastic Noise

Theorem 6.2 Let u be a ω-derivative mask with convergence order ρ and u l the corresponding ω-derivative l-skipping mask. Suppose that f : R -→ R is a C k function, with k ≥ 1 + ω, and

f (k) is bounded, α ∈]0, 1], K ∈ R * + and h ∈ R * + .
Let Γ be a digitization of f with step h and a stochastic noise ϵ h with expected value 0 and standard deviation σ(h) ≤ Kh α . Then for

k 0 = M in{k, 1+ρ} and l(h) = ⌊ h 1-α/k 0 ⌋ , the random variable ( 1 l(h) ω ∆ u l(h) ⋆ Γ ) (n)-f (ω) (nh)
has an expected value and a standard deviation which are O(h α-ωα/k 0 ).

The proof is similar to that of Theorem 6.1.

Examples of 1-derivative masks

It should be noted that all the masks defined in this section, which are digital 1-derivative (ω = 1) skipping masks, can be proven to have a multigrid convergence property either from Theorem 6.1 or from Theorem 6.2 (depending on the model of errors on values). We remind the reader that a more precise idea of the speed of convergence for the different masks can be found in Appendix A.2.

Binomial 1-Derivative Mask

Definition 7.1 (Binomial skipping masks) Let m and l be positive integers and i be an integer. The binomial 1-derivative l-skipping mask with diameter 2ml and cardinality 2m + 1 for discrete first order derivation is the finite sequence of rational numbers B

(1) l,2m-1 defined by B

(1)

l,2m-1 (i) = 1 l2 2m-1 ( ( 2m-1 m-1-i l ) - ( 2m-1 m-i l ) )
if l divides i and -ml ≤ i ≤ ml, and by B

(1) l,2m-1 (i) = 0 in any other cases.

Remark 7.1 Note that B (1) l,2m-1 is the convolution product ∆ + l ⋆B (0)
l,2m-1 . From Proposition 3.1, this is a general method for building a derivative mask from a smoothing mask.

Proposition 7.1 The masks B

(1) l,2m-1 are discrete 1-derivative masks of convergence orders 2.

Proof. We prove here the following equalities:

∑ i∈Z B (1) l,2m-1 (i) = 0 ∑ i∈Z iB (1) l,2m-1 (i) = -1 ∑ i∈Z i 2 B (1) l,2m-1 (i) = 0 ∑ i∈Z i 3 B (1) l,2m-1 (i) = 1-3m 2 l 2
Let us notice that for all k, we have:

∑ i∈Z i k B (1) l,2m-1 = l k-1 2 2m-1 ∑ j=m j=-m j k ( ( 2m-1 m-1-j ) - ( 2m-1 m-j
) )

and

∑ j=-1 j=-m j k ( ( 2m-1 m-1-j ) - ( 2m-1 m-j ) ) = ∑ j=m j=1 (-j) k ( ( 2m-1 m-1+j ) - ( 2m-1 m+j ) ) ∑ j=-1 j=-m j k ( ( 2m-1 m-1-j ) - ( 2m-1 m-j ) ) = ∑ j=m j=1 j k ( ( 2m-1 m-j ) - ( 2m-1 m-1-j ) )
Hence, for the even values of k, we have

∑ i∈Z i k B l,2m-1 (i) = 0.
For the odd values of k, we have:

∑ i∈Z i k B (1) l,2m-1 = p k-1 2 2m-1 ( ∑ j=m-1 j=-m j k ( 2m-1 m-1-j ) - ∑ j=m j=-m+1 j k ( 2m-1 m-j ) ) ∑ i∈Z i k B (1) l,2m-1 = p k-1 2 2m-1 ( ∑ j=m j=-m+1 (j -1) k ( 2m-1 m-j ) - ∑ j=m-1 j=-m j k ( 2m-1 m-j
) )

∑ i∈Z i k B (1) l,2m-1 = p k-1 2 2m-1 ∑ j=m j=-m ( (j -1) k -j k ) ( 2m-1 m-j ) For k = 1, we get ∑ i∈Z iB (1) l,2m-1 = -1.
For k = 3, we use previous results:

∑ i∈Z i k B (1) l,2m-1 = p k-1 2 2m-1 ∑ j=m j=-m (-3j 2 + 3j -1) ( 2m-1 m-j ) ∑ i∈Z i k B (1) l,2m-1 = p k-1 2 2m-1 (-3m.2 2m-2 + 3.2 2m-2 -2 2m-2 ) 2 Proposition 7.2 The masks B (1)
l,2m are discrete 1-derivative masks of convergence orders 1.

The proof could be derived from Remark 7.1 and Proposition 3.1, or we could follow the first steps of the proof of Proposition 7.1.

Taylor-Optimal 1-Derivative Mask

Definition 7.2 (Taylor-Optimal l-skipping masks) Let m and l be positive integers and i be an integer. The Taylor-optimal 1-derivative l-skipping mask with diameter 2ml and cardinality 2m + 1 for discrete first order derivation is the finite sequence of rational numbers

O (1) l,m defined by O (1) l,m (i) = (-1) i l i ( m | i l | ) ( m+| i l | m ) = (-1) i l i ( 2m m+| i l | ) ( 2m m )
if l divides i and -ml ≤ i ≤ ml, and by

O (1)
l,m (i) = 0 in any other cases.

Let us note that the complexity of the computation depends mainly on the parameter m and a few on the parameter l.

Proposition 7.3 O (1)
l,m are discrete 1-derivative mask with convergence order 2m .

The proof follows easily from Proposition 5.2.

B-spline 1-Derivative Mask

l,r (i) = 1 l ( N 0,r-1 ( r 2 + i l ) -N 0,r-1 ( r 2 + i l -1 )) if l divides i and -⌊ lr 2 + 1⌋ ≤ i ≤ ⌊ lr 2 + 1⌋ and S (1) l,r (i) = 0 otherwise.
This mask is called the B-spline uniform l-skipping 1-derivative mask with degree r -1.

Remark 7.2 Note that S

(1)

l,r = ∆ - l * S (0) l,r-1 , where S (0)
l,r-1 denotes the l-skipping smoothing mask associated to the uniform B-spline smoothing mask S (0) r-1 with degree r -1.

Proposition 7.4 The mask S

(1) l,r is a 1-derivative mask with convergence order 1.

The proof could be derived from Remark 7.2 and Proposition 3.1, but we present a hand made proof. Proof. The definition of a 1-derivative mask with convergence order 1 consists of two sum equalities. Hence the proof consists of two steps:

First Step. l ∑ i∈Z S (1) l,r (i) = ∑ i∈lZ N 0,r-1 ( r 2 + i l + 1 ) -N 0,r-1 ( r 2 + i l -1 ) = 0.
The last equality follows from change of index i ↔ (i -l) in one of the sums.

Second Step. l ∑ i∈Z i.S (1) l,r (i) = ∑ i∈lZ i. [ N 0,r-1 ( r 2 + i l ) -N 0,r-1 ( r 2 + i l -1 )] = ∑ i∈lZ i.N 0,r-1 ( r-1 2 + i l ) - ∑ i∈lZ (i + l).N 0,r-1 ( r-1 2 + i l ) = -l ∑ i∈lZ N 0,r-1 ( r-1 2 + i l ) = -l 2
The following methods for building derivative masks from a smoothing mask u of order ρ are worth noting:

1. Let v be the sequence defined by v(0)

= 1 ∑ j∈Z,j̸ =0 u(j) ∑ j∈Z,j̸ =0 1 j u(j) and v(i) = -u(i) i ∑ j∈Z,j̸ =0 u(j)
for i ̸ = 0. Then, v is a 1-derivative mask with order 1 + ρ if and only if u(0) = 0.

2. From Proposition 3.1 and Definition 3.2, we can easily define and (ω + 1)-derivative mask from any ω-derivative mask by covolution with ∆ + l .. Evaluating the errors needs to evaluate

∑ i∈Z |i k B (1) l,2m-1 (i)| and ∑ i∈Z |i k O (1)
l,m (i)| (see the proof of Theorem 6.1). See Appendix A.4 for the most useful such estimations.

8 The cases of 2-Derivatives and Parametric Curves

Exemples of 2-derivative masks

Concerning second order derivatives, we can consider similar l-skipping masks. It is worth noting the following general method: Let u and v be two 1-derivative masks, then u ⋆ v is a 2-derivative mask (see Proposition 3.1). Definition 8.1 Let m and l be positive integers and i an integer. The Taylor-optimal l-skipping mask with diameter 2ml and cardinality 2m + 1 for discrete second order derivation is the finite sequence of rational numbers

O (2) l,m defined by O (2) l,m (i) = -2 li O (1) l,m (i) if l divides i and -ml ≤ i ≤ ml and i ̸ = 0 and O (2) l,m (0) is such that ml ∑ i=-ml O (2)
l,m (i) = 0, and 0 for all other values of i.

Lemma 8.1 O (2)
l,m is a l-skipping 2-derivative mask with convergence order 1 + 2m. The proof follows easily from Proposition 5.2. Definition 8.2 Let m be a positive integers and i be an integer.

The binomial smoothing mask with diameter (2m+2)l and cardinality 2m+3 is the finite sequence of rational numbers

( B (2) l,2m-1 )
defined by B

(2)

l,2m-1 (i) = 1 l 2 2 2m-1 ( ( 2m-1 m+1-i l ) -2 ( 2m-1 m-i l ) + ( 2m-1 m-1-i l ) )
if -(m + 1)l ≤ i ≤ (m + 1)l and l divides i, and by B

(2) l,2m-1 (i)) = 0 in any other cases. See Appendix A.5 for estimations of constants in the convergence order.

Parametrized Curves

Let k ≥ 1. We assume that a planar closed C k -parametrized curve C (i.e. the parametrization is periodic) is given together with a family of parametrized digital curves (Σ h ) with Σ h contained in a tube with radius H(h) around C. A parametrization of C is denoted by g = (g 1 , g 2 ) : R → R 2 , where g 1 and g 2 are C 1 and periodic with a common period. The parametrization of Σ h is denoted by

Γ h = (Γ 1 , Γ 2 ) : Z → Z 2 .

Binomial smoothing

Let P = (P 0 , . . . , P n ) be a (n + 1)-tuple of control points. Let us denote a Bézier curve

B P (t) = ∑ n i=0 ( n i ) t i (1 -t) n-i .P i ,
where t ∈ [0, 1] (See the Section 9.2). Thus we have

B (Γ(i 0 -m+ζ),...Γ(i 0 +m)) ( 1 2 ) = (( ∆ 0 B (0) 2m-ζ (Γ 1 ) ) (i 0 ), ( ∆ 0
B (0) 2m-ζ (Γ 2 ) ) (i 0 )
)

where ζ = 0 or 1.

Discrete tangent

Here we estimate the tangent at a point of C by a digital tangent at a not too far point of Σ h .

To this purpose, we define the digital derivative in a natural way with the digital derivative of the coordinates. Note that this definition needs no hypothesis concerning the parametrized curve. Moreover, we could theoretically use different masks for the two coordinates. For sake of simplicity, we shall use the same for both.

Definition 8.3 Let u be a 1-derivative mask. Let Γ = (Γ 1 , Γ 2 ) : Z → Z × Z be a digital parametrization of a digital curve Σ.

The digital derivative of Γ at point Γ(n) is ((

∆ 1 u (Γ 1 ) ) (n), ( ∆ 1 u (Γ 2 ) ) (n) )
It is denoted by (∆ 1 u (Γ)) (n). When this vector is nonzero, the corresponding discrete tangent at Γ(n) is the real line going through Γ(n) directed by (∆ 1 u (Γ)) (n).

Let us notice that the asymmetric l-skipping binomial case corresponds to the derivative of the Bézier curve with control points (Γ(i 0 -(m -1)l), ..., Γ(i 0 + ml)) computed at the medialpoint of parameter 1 2 .

Tangent Estimation

Theorem 8.1 Let u be a smoothing mask with convergence order ρ. Let g = (g 1 , g 2 ) be a C k parametrization of a closed curve Σ.

Let k 0 = inf {k, 1 + ρ}, α ∈]0, 1], K ∈ R * + and h ∈ R * + with h ≤ 1 2 . Let Γ h = (Γ h,1 , Γ h,2 ) : Z → Z × Z be a digital parametrization of Σ. Suppose that for all i we have ∥g(ih) -h.Γ h (i)∥ ∞ ≤ Kh α . Suppose that g (k 0 ) c is bounded for c = 1, 2. Then ∥∆ 1 u (Γ h )(n) -g ′ (nh)∥ ∞ is less than h k 0 -1 M ax { ∥g (k 0 ) 1 ∥ ∞ , ∥g (k 0 ) 2 ∥ ∞ } k 0 ! ( ∑ i∈Z |i k 0 u(i)| ) + K h 1-α ( ∑ i∈Z |u(i)| )
The proof follows directly from Lemma 4.2 and Lemma 4.1.

Theorem 8.2 Let u be a 1-derivative mask with convergence order ρ and u l the corresponding 1-derivative l-skipping mask. Under the hypothesis of Theorem 8.1, for k ≥ 1 and l(h

) = ⌊ h -1+α/k 0 ⌋ we have ∆ 1 u l(h) (Γ h )(n) -g ′ (nh) ∞ ∈ O(h α-α/k 0 )
, and for k = 1 and

l(h) = ⌊ h -1+α/2 ⌋ , we have ∆ 1 u l(h) (Γ h )(n) -g ′ (nh) ∞ ∈ o(1)
The proof is straightforward. The hypothesis of Theorem 8.1 is stronger than the digitization being simply contained in a tube, because not only the discrete curve must be close to the continuous curve, but the parametrization of the digitization must be close to the parametrization of the continuous curve. This was the reason for introducing Pixel-Length parametrizations in [START_REF] Esbelin | Convergence of binomial-based derivative estimation for c 2 noisy discretized curves[END_REF].

Discrete curvature

We estimate the curvature at a point of C by a digital curvature at a not too far point of Σ h . The definition of the digital curvature coincides with the classical formula for continuous curves.

Definition 8.4 Let u and v be respectively a 1-derivative mask and 2-derivative mask. Let Σ be a digital parametrized curve of parametrization Γ = (Γ 1 , Γ 2 ). The corresponding discrete curvature at Γ(n) is

(∆ 2 v (Γ 1 )) (n) (∆ 1 u (Γ 2 )) (n) -(∆ 1 u (Γ 1 )) (n) (∆ 2 v (Γ 2 )) (n) ( ((∆ 1 u (Γ 1 )) (n)) 2 + ((∆ 1 u (Γ 2 )) (n))
2 ) 3/2 when the vector ∆ 1 u (Γ)(n) is non-zero. It will be denoted by C u,v (Γ)(n). Notice that this definition needs no more hypothesis concerning the parametrized curve.

Curvature Estimation

We shall denote γ g (t) the curvature of a curve parametrized by g at the point of parameter t. The following convergence theorem required some additional hypothesis concerning the regularity of the considered point.

Theorem 8.3 Let u be a 1-derivative mask with convergence order ρ(u) and u l the corresponding 1-derivative l-skipping mask. Let v be a 2-derivative mask with convergence order ρ(v) and v l the corresponding 2-derivative l-skipping mask. Let g = (g 1 , g 2 ) be a C k parametrization of a simple closed curve, with k ≥ 2.

Let Γ = (Γ 1 , Γ 2 ) : Z → Z × Z be a digital parametrization of a digital curve Σ. Let k 0 = inf {k, 1 + ρ(u), 1 + ρ(v)}, α ∈]0, 1], K ∈ R *
+ and h ∈ R * + and h ≤ 1 2 . Suppose that for all i we have ∥g(nh) -h.Γ(n)∥ ∞ ≤ Kh α . Suppose that g (k 0 ) c is bounded for c = 1, 2. Suppose that for t ∈ [t 1 , t 2 ], the vector g ′ (t) ̸ = 0. Then, there is some h 0 > 0 and two functions n 1 , n 2 : (0; +∞) → Z such that for l

(h) = ⌊ h -1+α/k 0 ⌋ , for h ≥ h 0 and n ∈ [n 1 (h), n 2 (h)] ∩ Z, we have C u l(h) ,v l(h) (Γ)(n) -γ g (nh) ∈ O(h α-2α/k 0 )
Definition 9.2 The uniform l-scaled B-spline function with degree r associated to Γ is defined by: ( S (0) l,r (Γ)

) (t) = ∑ s∈Z Γ(ls)N s,r,l (t)
Note that this last function can be interpreted a regular uniform B-spline functions control points (Γ(ls)) s∈Z with difference between successive nodal points equal to l. Definition 9.3 Let Φ : Z -→ E be a sequence, or Φ : R -→ E be a function. We denote by τ (Φ) the right shift of Φ which to s associates (τ (Φ)) (s) = Φ(s + 1).

The following property follows immediately from the shift property in Lemma 5.1.

Property 9.1 ( S (0)
l,r (τ l (Γ))

) (t) = ( τ l (S (0)
l,r (Γ))

) (t)

Derivative of the B-splines Curves with Skipped Control Points

As far as the Cox de-Boor functions are concerned, we have the classical expression for the derivative:

N ′ s,r (t) = (N s,r-1 -N s+1,r-1 ) From the definition of N s,r,l (t) we get N ′ s,r,l (t) = 1 l N ′ s,r ( t l ) = 1 l (N s,r-1,l (t) -N s+1,r-1,l (t)) (6)
Proposition 9.1 For r ∈ N * we have:

d dt ( S (0) l,r (Γ) ) (t) = 1 l ( S (0) l,r-1 (Γ) -S (0) l,r-1 (τ -l (Γ)) ) Proof. d dt ( S (0) l,r (Γ) ) (t) = ∑ s∈Z Γ(ls)N ′ s,r,l (t) = ∑ s∈Z Γ(ls) 1 l (N s,r-1,l (t) -N s+1,r-1,l (t)) = ∑ s∈Z Γ(ls) 1 l (N s,r-1,l (t)) - ∑ s∈Z Γ(l(s -1)) 1 l (N s,r-1,l (t)) = 1 l ∑ s∈Z ( Γ(ls) -τ -l Γ(ls) ) (N s,r-1,l (t)) 2
Definition 9.4 Let Φ : Z -→ E be a sequence, or Φ : R -→ E be a function. We denote

( ∆ - l (Φ) ) (s) = 1 l (Φ(s) -Φ(s -l)).
Notation 9.1 For ω ≤ r -1, we denote by S (ω) l,r (Γ) the function defined on R by: ( S (ω) l,r (Γ)

) (t) = d ω dt ω ( S (0) l,r (Γ) ) (t)
Notation 9.4 For ω ∈ {0, . . . , r -1}, we denote by B

(ω) l,r (Γ) the mask on R defined by

B (ω) l,r (n) = ( ∆ + l ) ω (B l,r-ω ) (n)
Note that this definition is coherent with Definition 7.1 and that B

(ω) l,r is an ω-derivative mask.

Theorem 9.2 For r ≥ 1 and ω ≤ r -1, we have

( B (ω) l,r (Γ) ) ( l(n + 1 2 ) ) = r! (r -ω)! ( B (ω) l,r * Γ ) ( l(n -⌊ r 2 ⌋)
)

Proof. The result follows by induction from Proposition 9.6 and Proposition 9.7. 2 10 Experiments 10.1 Smoothing Figure 1 shows an example of smoothing a sampled function, with different models of noise of the sample values. Figure 2 shows the maximal error for smoothing with uniform noise, for different values of the resolution 1 h . We see that the Taylor-Optimal masks, which are optimal in the absence of noise, perform poorly for removing noise. The result of "smoothing" in this case is even more noisy than the original.

The best practical kernels seem to be that Binomial kernels. Although Proposition 5.1 indicates that the convergence order of the binomial masks with even or odd diameter is different, the practical performance for smoothing appears similar.

Derivatives

Figure 3 shows that binomial 1-derivative masks have a better practical error than the Tayloroptimal one in various cases with noise, which confirms the theoretical bounds. The results are for a uniform noise with maximum absolute value 1 4 h 0.8 (i.e. α = 0.8). Mask cardinality r = 11. Skipping step for l =⌋h -0.5 ⌊ and for l =⌋h -0.6 ⌊ Figure 4, also confirms the result that the Taylor-Optimal 1-derivative masks perform better in the absence of noise, but similar, or even worse in the presence of very moderate noise (close to round-off quantisation errors). We see that the performance of the Binomial masks have a closer dependency to the resolution, but still have a broad spectrum of good performance.

We also see on both Figure 3 and Figure 4 that a good estimation can be obtained by a mask using only 11 (possibly noisy) values. This shows the strength, as far as computational complexity is concerned, of the method with skipping steps compared to using contiguous samples. Figure 5 shows that, also with only 40 points, the second derivative can be estimated with good precision in the round-off quantization value errors case. The figure also shows the influence of the skipping step. Bellow some threshold power of h, the method just does not converge.

Conclusion

This work is a milestone in our work in progress to define a new paradigm for numerical analysis and calculus, to avoid the shortcomings of floating point computations. It can be easily extended to any dimension by tensor product, which leads to convergence results by application of the one-dimensional case. The interpretation of digital convolutions using smooth curves can be generalized to many convolution kernels, including the Gaussian kernels, for which the smooth curve can be chosen C ∞ . Even though we are primarily interested in integer-only or exact computations, this is an interesting limit case, which can be of interest to many scientists using more traditional approaches to numerical analysis ans functional analysis.

We will build on the general notion of a digital derivative to define other tools for digital analysis, such as integrals, partial differential equations, convex optimization, harmonic analysis, and non linear optimization. Some of this work is already in progress as we write these lines.

We also plan to develop a dedicated library, for efficient integer-only computations. The library should be generic and support many formats for numbers, including long-integers and floating points. More experiments on n-dimensional derivative and curvature estimation will be performed in this process. Even-though long-integers come with an overhead, we believe that exact calculations will prevail over cheap approximations in the long run.

A Choosing a Mask

A.1 Criteria for Choosing a Mask

Having read all of the paper until the conclusion, the reader may still be wondering: "Why bother with so many masks?". In fact, there is no single bullet rule for choosing a mask. For now, any user of our methods has to assess what we call "business dependant" or "application dependant" hypothesis, about available resolutions and noise model.

If you have low resolution and low noise (a few pixels), you should probably use the Taylor Optimal masks (see Figure 4). If you have a fair resolution and significant noise, you will probably be better off with a Binomial mask.

The B-spline masks have not yet been completely assessed, and they are here for their potential theoretical interest, which deserves the proofs of multigrid convergence to be included.

In the next section, we attempt to provide finer criteria to choose masks for your favourite application.

A.2 Moments of The Main Considered Masks

In Lemma 4.1, an upper bound is given for the sampling error |ES ω (f, h, Γ, u, n)|, as a function of the infinite norm of f (k 0 ) , the digitization step, and a constant which depends on u: ∑ i∈Z |i k 0 u(i)|. For k 0 = 1 + ρ, where ρ is the convergence order of the digital ω-derivative mask u, we call this last constant, the moment of the mask.

This constant gives an indication on the speed of convergence (main asymptotic term) of the sampling error the digital ω-derivative estimation given by u for C 1+ρ functions. We can also consider this constant for k 0 ≤ ρ to assess the quality of estimation on functions which are not C 1+ρ .

Similarly, in Lemma 4.2 and in the proof of Theorem 6.1, the main asymptotic term in the error on values |EV ω (f, h, Γ, u, n)| depends on a constant ∑ i∈Z |u(i)|. We also asses this values for some masks.

As we advocate the use of masks of small support, we give here schedules for the coefficients values.

A.3 Smoothing Masks

  (i) = 0 in any other cases. The binomial smoothing mask with diameter 2m and cardinality 2m+1 is the finite sequence of dyadic rational numbers B

  ) = 0 in any other cases.

  Proposition 5.1 The masks B (0) r

  and by O (0) 2m (i) = 0 in any other cases. Proposition 5.2 The mask O

  Definition 7.3 (B-spline Derivative Kernel) Let us consider two integers r ≥ 0 and l > 0. We introduce a convolution mask S (1) l,r , with support in {-⌊ lr 2 + 1⌋, . . . , ⌊ lr 2 + 1⌋}, based on the uniform Cox-de Boor functions with degree r -1 on R. The mask values are defined by S (1)

Figure 1 :

 1 Figure 1: Smoothing operators applied to the sin function, digitization step h = 1 50

Figure 2 :

 2 Figure 2: Smoothing operators applied to the sine function, uniform noise of maximum absolute value 2h.

Figure 3 :

 3 Figure 3: Maximal error 1-derivative operators applied to the square function, for different values of the resolution 1 h from 200 to 1000. The results are for a uniform noise with maximum absolute value 1 4 h 0.8 (i.e. α = 0.8). Mask cardinality r = 11. Skipping step for l =⌋h -0.5 ⌊ and for l =⌋h -0.6 ⌊

Figure 4 :

 4 Figure 4: Maximal error for 1-derivative operators applied to the sin function, uniform noise of maximum absolute value 2h, mask cardinality r = 11 and l = 100

Figure 5 :

 5 Figure 5: Maximal error of 2-derivative operators with binomial skipping masks, on the sin function, with only round-off quantization value errors. All convolution values are computed with 40 samples, for different values l = h 0.7 , l = 40h 1/3 , l = h 0.8 , l = 15h 1/3 , l = h 0.9

Figure 6 :

 6 Figure 6: Digital curvature around a noisy digitization of an ellipse, uniform noise of maximum absolute value 3h, h = 1 500 , Taylor-optimal mask, various values for l

  2m (i)| 0.157 0.089 0.051 0.029 0.017 0.0.010 0

  Aim of the work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Errors Related to Sampling and to Input Values . . . . . . . . . . . . . . . . . . 8 4.2 Upper Bound for the Sampling Error . . . . . . . . . . . . . . . . . . . . . . . . 8 4.3 Upper Bound for the Values Error . . . . . . . . . . . . . . . . . . . . . . . . . . 9 Binomial 1-Derivative Mask . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 7.2 Taylor-Optimal 1-Derivative Mask . . . . . . . . . . . . . . . . . . . . . . . . . 15 7.3 B-spline 1-Derivative Mask . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 Continuous B-Splines Revisited . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 9.1.1 Scaled B-Splines Curves Associated to a Sequence . . . . . . . . . . . . 19 9.1.2 Derivative of the B-splines Curves with Skipped Control Points . . . . Smoothing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 10.2 Derivatives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 10.3 Curvature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 Criteria for Choosing a Mask . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 A.2 Moments of The Main Considered Masks . . . . . . . . . . . . . . . . . . . . . . 29 A.3 Smoothing Masks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 A.4 First Order Derivative Masks . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 A.5 Second Order Derivatives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

	9 Piecewise Polynomial Functions and Convolutions	19
	9.1 . 20
	9.1.3 Relationship between scaled B-spline and convolutions . . . . . . . . . . 21
	9.2 Bézier Curves with Skipped Control Points . . . . . . . . . . . . . . . . . . . . . 21
	9.2.1 Bézier Curves and Bernstein polynomials Basics . . . . . . . . . . . . . . 21
	9.2.2 Scaled Bézier Function Associated to a Sequence . . . . . . . . . . . . . . 22
	9.2.3 Derivative of the Scaled Bézier Function . . . . . . . . . . . . . . . . . . 23
	9.2.4 Binomial Convolutions and Scaled Piecewise Bézier . . . . . . . . . . . . 24
	10 Experiments	25
	1 Introduction 1.1 2 Digitization, Quantization, Noise Models 10.1 11 Conclusion	1 5 28
	A Choosing a Mask	28
	3 Digital Derivatives for Functions 4 Basic Error Decomposition and Upper Bounds 4.1 5 Smoothing masks 5.1 Binomial Smoothing Masks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 6 8 10 A.1 2 Digitization, Quantization, Noise Models
	5.2 Taylor-Optimal Smoothing Masks . . . . . . . . . . . . . . . . . . . . . . . . . . 11
	5.3 B-splines Smoothing Masks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
	6 Skipping Masks: Cheap Multigrid Convergence	13
	6.1 Uniform Multigrid Convergence with Uniform Noise or Bias . . . . . . . . . . . 13
	6.2 Stochastic Multigrid Convergence with Stochastic Noise . . . . . . . . . . . . . . 14
	7 Examples of 1-derivative masks	14
	7.1 8 The cases of 2-Derivatives and Parametric Curves	16
	8.1 Exemples of 2-derivative masks . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
	8.2 Parametrized Curves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
	8.2.1 Binomial smoothing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
	8.2.2 Discrete tangent . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
	8.2.3 Tangent Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
	8.2.4 Discrete curvature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
	8.2.5 Curvature Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

Table 1 :

 1 Moments of the mask O

	(0)
	2m

Table 2 :

 2 Moments of the mask B

	(0)
	2m-1

Table 3 :

 3 Moments of the mask SIn practice, m will be 3 to 7. Hence we first give a table of values for such m. 2 7.32l 2 10.2l 2 13.3l 2 16.8l 2 20.4 l 2 24.3l 2

	(0)
	r

Table 4 :

 4 ∑ i∈Z |i k O

	(1) l,m (i)|

Table 5 :

 5 ∑ i∈Z |i k B

	(1) l,2m-1 (i)|

1 This work is supported by Agence Nationale pour la Recherche (A.N.R.), KIDICO (ANR 2010 BLANC 0205 02) project

Proof. Let us introduce the function defined on R 2 -{(0, 0)} × R 2 by

Hence C u l(h) ,v l(h) (Γ)(n) = F (∆ 1 u l(h) (Γ)(n), ∆ 2 v l(h) (Γ)(n)) and γ g (nh) = F (g ′ (nh), g"(nh)).

The differential of F is easily seen to be bounded on any compact K in R 2 -{(0, 0)} × R 2 . Let M K be one of its upper bounds on K. Then for a and b in K, we have ∥F (a) -

We claim that, under the hypothesis of the theorem, we know that there is a compact K 1 in R 2 -{(0, 0)} such that for t ∈ [t 1 , t 2 ] ⊂ R, the derivative g ′ (t) lies in K 1 . Indeed K 1 = {x ∈ R 2 ; α ≤ ∥x∥ ≤ β}, where α = Inf {∥g ′ (t)∥; t ∈ [t 1 ; t 2 ]} and β = M ax{∥g ′ (t)∥; t ∈ [t 1 ; t 2 ]} is convenient.

Then from Theorem 8.2, we can find a compact K 2 in R 2 -{(0, 0)} and a fixed h 1 such that for h ≥ h 1 and n ∈ [⌈ t 1 h ⌉, ⌊ t 1 h ⌋] ∩ Z the derivatives g ′ (nh) and (∆ 1 u l(h) (Γ)(n) ly in K 2 . Moreover, under previous conditions, the second order derivatives g"(nh) and

From Theorem 6.1, there is a constant χ and a fixed h 2 such h ≥ h 2 , we have

where the N s,r function is the Cox De-Boor function (see Definition 5.3).

Definition 9.1 Let l ∈ N * + and r ∈ N. First, for s ∈ Z we define the l-scaled B-spline function N s,r,l by:

With these notations, Proposition 9.1 can be reformulated as

))

Hence by induction, we get:

Proposition 9.2 For ω ≤ r -1, we have

))

9.1.3 Relationship between scaled B-spline and convolutions Proposition 9.3 For n ∈ Z, we have

) ( l

Notation 9.2 For ω ∈ {0, . . . , r -1}, we denote

Note that, from Proposition 3.1, the mask S (ω) l,r is an ω-derivative mask.

Theorem 9.1 For n ∈ Z and ω ∈ {0, . . . , r -1}, we have

) ( l

9.2 Bézier Curves with Skipped Control Points

Bézier Curves and Bernstein polynomials Basics

Definition 9.5 Let r ∈ N * . For i ∈ 0, . . . , r, we consider the R-valued polynomial with degree r defined by

These well-known polynomials are called the Bernstein polynomials with degree r.

In the sequel, unless otherwise specified, we shall say Bernstein polynomials or Bernstein functions as a shorthand for Bernstein polynomials. The Bernstein polynomials with degree r constitute a basis of the vector space of polynomials with degree less than or equal to r.

From the Pascal formula for binomial coefficients, we derive a similar formula about Bernstein polynomials:

Definition 9.6 Let E be an R-module and let P = (P 0 , . . . , P r-1 ) be a sequence of point in E. Let B P : 1 -→ E be the Bézier Curve with control points P. Then we have:

We can compute the derivative of a Bézier Curve as follows. Let B P : [0, 1] -→ E be the Bézier Curve with control points P 0 , . . . , P r-1 . Then we have:

which is the Bézier curve with control points P ′ 1 , . . . , P ′ n-1 , where P ′ i = (r -1)(P i+1 -P i ).

Scaled Bézier Function Associated to a Sequence

Definition 9.7 Let l ∈ N * + . We introduce the s-shifted l-scaled Bernstein polynomials by:

In other words, the value B i,s,r,l (t) can be non-zero only for s = ⌊ t l ⌋. In the remainder of this section, E denotes an R-module and (Γ(s)) s∈Z is a sequence with values in E. Definition 9.8 Let l ∈ N * + . The l-scaled (piecewise) Bézier curve with degree r-1 associated to Γ is defined over R by:

Note that in the previous definition, due to the definition of B i,s,r,l , for a given value of i, only one value of s (namely s = ⌊ t l ⌋ + r -i) contributes to the double sum

, so that, in fact, at most r + 1 terms are non-zero.

)

The last equality follows from Equation 7. Taking into account that the only value of s for which t l ∈ [s, s + 1[, and then by changing the index i to i + 1 in the sum, we get:

)

The indices i = -1 and i = r yielding a zero term because out of range for the Bernstein polynomials. 2

Derivative of the Scaled Bézier Function

As far as Bernstein polynomials are concerned, we have the classical formula:

We derive from this that, for s = ⌊t⌋, we have

so that Proposition 9.5 For r ≥ 1, the curve B (0) l,r (Γ) is C r-1 on R and we have:

Proof. First we prove the result for all t ∈ R\Z, on which the curve B (0) l,r (Γ) is easily seen to be polynomial.

) Now, as we can see that B (0) l,1 (Γ) is C 0 on R, so that the result follows by induction on r ≥ 2. 2 Definition 9.9 Let Φ : Z -→ E be a sequence, or Φ : R -→ E be a function. We denote

Notation 9.3 For ω ∈ {0, . . . , r -1}, we denote by B

(ω) l,r (Γ) the function on R defined by

Therefore, Proposition 9.5 can be restated as:

l,r (Γ)

The following immediately follows by induction:

l,r (Γ)

Binomial Convolutions and Scaled Piecewise Bézier

Proposition 9.7

We distinguish two cases depending on the parity of r: First Case: If r is even. Say r = 2m with m ∈ N. We have:

In order to compare both masks, let us consider a C 2 function f and m = 5. In the binomial case, the error bound is 7

l(h) . In the optimal case, the error bound is 11

l(h) which is worse. For a C 3 function, the error bound remains the same in the binomial case, but in the optimal case, it become 25 24 ||f

l(h) . The binomial mask seems to be competitive in the generic case. Hence, we give here general formulas for this case:

Now for even k, we have

) ) .

For k = 0, we have

) .

For k ≥ 2, we have

)

For k = 2, we have

∑ j=m j=1 (2j -1)

) 2

A.5 Second Order Derivatives


Michaël Gutnic

Ioana Paun

Eric Sonnendrücker

Vlasov simulations on an adaptive phase-space grid

come

Vlasov simulations on an adaptive phase-space grid

M. Gutnic 1 , I. Paun 1 , E. Sonnendrücker 1 1

IRMA, Université Louis Pasteur, Strasbourg, France

The numerical resolution of the Vlasov equation is usually performed by particle methods (PIC) which consist in approximating the plasma by a finite number of particles. The trajectories of these particles are computed from the characteristic curves given by the Vlasov equation, whereas selfconsistent fields are computed on a mesh of the physical space. This method allows to obtain satisfying results with a relatively small number of particles. However, it is well known that, in some cases, the numerical noise inherent to the particle method becomes too important to have an accurate description of the distribution function in phase space. To remedy to this problem, methods discretizing the Vlasov equation on a mesh of phase space have been proposed [START_REF] Cheng | The integration of the Vlasov equation in configuration space[END_REF][START_REF] Sonnendrücker | Simulation of heavy ion beams with a semi-Lagrangian Vlasov solver[END_REF][START_REF] Sonnendrücker | The Semi-Lagrangian Method for the Numerical Resolution of Vlasov Equations[END_REF].

The major drawback of Vlasov methods using a uniform and fixed mesh is that their numerical cost is high, which makes them rather inefficient when the dimension of phase-space grows. For this reason we are investigating a method using an adaptive mesh. The adaptive method is overlayed to a classical semi-Lagrangian method which is based on the conservation of the distribution function along particle trajectories. The phase-space grid is updated using a multiresolution technique.

The model we consider throughout this paper is the nonrelativistic Vlasov equation coupled self-consistently with Poisson's equation. It reads

∂f ∂t + v • ∇ x f + q m (E + v × B) • ∇ v f = 0, (1)
the self electric field E is computed from Poisson's equations

-0 ∇ 2 φ = ρ(x, t) = q f (x, v, t) dv, E = -∇φ.
The magnetic field is external and considered to be known.

In the present work, we have chosen to introduce a phase-space mesh which can be refined or derefined adaptively in time. For this purpose, we use a technique based on multiresolution analysis which is in the same spirit as the methods developed in particular by S. Bertoluzza [START_REF] Bertoluzza | An adaptive collocation method based on interpolating wavelets[END_REF], A. Cohen et al. [START_REF] Cohen | Fully adaptive multiresolution finite volume schemes for conservation laws[END_REF] and M. Griebel and F. Koster [START_REF] Griebel | Adaptive wavelet solvers for the unsteady incompressible Navier-Stokes equations[END_REF]. We represent the distribution function on a wavelet basis at different scales. We can then compress it by eliminating coefficients which are small and accordingly remove the associated mesh points. Another specific feature of our method is that we use an advection in physical and velocity space forward in time to predict the useful grid points for the next time step, rather than restrict ourselves to the neighboring points. This enables us to use a much larger time step, as in the semi-Lagrangian method the time step is not limited by a Courant condition. Once the new mesh is predicted, the semi-Lagrangian methodology is used to compute the new values of the distribution function at the predicted mesh points, using an interpolation based on the wavelet decomposition of the old distribution function. The mesh is then refined again by performing a wavelet transform, and eliminating the points associated to small coefficients.

This paper is organized as follows : we shall first present the context of multiresolution analysis for adaptive simulations, then describe the adaptive semi-Lagrangian method and validate our method on typical simulation in beam physics.

Multiresolution analysis

The semi-Lagrangian method consists mainly of two steps, an advection step and an interpolation step. The interpolation part is performed using for example a Lagrange interpolating polynomial on a uniform grid. Thus interpolating wavelets provide a natural way to extend this procedure to an adaptive grid in the way we shall now shortly describe.

For simplicity, we shall restrict our description to the 1D case of the whole real line. It is straightforward to extend it to periodic boundary conditions and it can also be extended to an interval with Dirichlet boundary conditions. The extension to higher dimension is performed using a tensor product of wavelets and will be addressed at the end of the section.

For any value of j ∈ Z, we consider a uniform grid G j of step 2 -j . The grid points are located at x j k = k2 -j . This defines an infinite sequence of grids that we denote by (G j ) j∈Z , and j will be called the level of the grid.

In order to map a function defined on the grid from one level to the next or the previous, we define a projection operator and a prediction operator. Consider two grid levels G j and G j+1 and discrete values (of a function) denoted by (c j k ) k∈Z and (c j+1 k ) k∈Z . Even though we use the same index k for the grid points in the two cases, there are of course twice as many points in any given interval on G j+1 as on G j . Using the terminology in [START_REF] Cohen | Fully adaptive multiresolution finite volume schemes for conservation laws[END_REF], we then define the projection operator

P j j+1 : G j+1 → G j , c j+1 2k → c j k
, which is merely a restriction operator, as well as the prediction operator

P j+1 j : G j → G j+1 , such that c j+1 2k = c j k , c j+1 2k+1 = P 2N -1 (x j+1 2k+1 ),
where P 2N -1 stands for the Lagrange interpolation polynomial of odd degree 2N -1 centered at the point (x j+1 2k+1 ). Using the just defined prediction operator, we can construct on G j a subspace of L 2 (R) that we shall denote by V j , a basis of which being given by (ϕ j k ) k∈Z such that ϕ j k (x j k ) = δ kk where δ kk is the Kronecker symbol. The value of ϕ j k at any point of the real line is then obtained by applying, possibly an infinite number of times, the prediction operator.

In the wavelets terminology the ϕ j k are called scaling functions. Let us denote by

ϕ = ϕ 0 0 . Notice that [-N + 1, N -1], ϕ j k (x) = ϕ(2 j x -k).
The scaling functions have the following noteworthy properties: 1) they have compact support, ϕ vanishes outside, 2) they are interpolating as by construction ϕ(0) = 1 and ϕ(k) = 0 if k = 0, 3) All polynomials of degree less than N -1 can be obtained exactly as linear combinations of the ϕ 0 k , 4) Change of scale: ϕ can be expressed using the half scale:

ϕ(x) = N -1 -N +1 h l ϕ(2x -l),
the h l are linked to the projection operator P j+1 j , as c j+1 k = h l-2k c j l . As V j ⊂ V j+1 , there exists a supplementary of V j in V j+1 that we shall call the detail space and denote by W j :

V j+1 = V j ⊕ W j .
The construction of W j can be made in the following way: an element of V j+1 is characterized by the sequence(c j+1 k ) k∈Z and by construction we have

c j k = c j+1 2k .
Thus, if we define

d j k = c j+1 2k+1 -P 2N +1 (x j+1 2k+1
), where P 2N +1 is the Lagrange interpolation polynomial by which the value of an element of V j at the point (x j+1 2k+1 ) can be computed, d j k represents exactly the difference between the value in V j+1 and the value predicted in V j . Finally, any element of V j+1 can be characterized by the two sequences

(c j k ) k of values in V j and (d j k ) k of details in W j .
Moreover this strategy for constructing W j is particularly interesting for adaptive refinement as d j k will be small at places where the prediction from V j is good and large elsewhere, which gives us a natural refinement criterion. Besides, there exists a function ψ, called wavelet such that {ψ j k = 2 j/2 ψ(2 j x -k)} k∈Z is a basis of W j . In practise, for adaptive refinement we set the coarsest level j 0 and the finest level j 1 , j 0 < j 1 , and we decompose the space corresponding to the finest level on all the levels in between:

V j 1 = V j 0 ⊕ W j 0 ⊕ W j 0 +1 ⊕ • • • ⊕ W j 1 -1 .
A function f ∈ V j 1 can then be decomposed as follows

f (x) = +∞ l=-∞ c j 0 l ϕ j 0 l (x) + j 1 -1 j=j 0 +∞ l=-∞ d j l ψ j l (x),
where the (c j 0 l ) l are the coefficients on the coarse mesh and the (d j l ) l the details at the different level in between.

2D generalisation:

In two dimensions, the prediction operator which is the base of the multiresolution analysis can be defined as a tensor product of 1D operators. In practice, one needs to consider three distinct cases. (see figure 1 for notations):

1. Refinement in x (corresponds to the points c j+1 2k 1 +1,2k 2 and c j+1 2k 1 +1,2k 2 +2 ): we use the 1D prediction operator in x for fixed k 2 .

Refinement in v (corresponds to the points c j+1

2k 1 ,2k 2 +1 and c j+1 2k 1 +2,2k 2 +1 ): we use the 1D prediction operator in v for fixed k 1 .

Refinement in x and v (corresponds to the point c j+1

2k 1 +1,2k 2 +1 ): we first use the 1D prediction operator in v for fixed k 1 to construct the necessary points to apply the 1D prediction operator in x for fixed k 2 that is the applied.

The corresponding wavelet bases are respectively of type ψ(x)ϕ(v), ϕ(x)ψ(v) et ψ(x)ψ(v) where ϕ and ψ are respectively the 1D scaling function and wavelet. Indeed if we have a decomposition like

V j+1 × V j+1 = V j × V j + V j × W j + W j × V j + W j × W j .
we obtain a 2D wavelet decomposition of the following type

f (x, v) = k 1 ,k 2 c j 0 k 1 ,k 2 ϕ j 0 k 1 (x) ϕ j 0 k 2 (v) + j 1 -1 j 0 d row,j k 1 ,k 2 ψ j k 1 (x) ϕ j k 2 (v) + d col,j k 1 ,k 2 ϕ j k 1 (x) ψ j k 2 (v) + d mid,j k 1 ,k 2 ψ j k 1 (x) ψ j k 2 (v) . ( 2
)
nD generalization : This will be useful for three, four, five or six dimensional phase space as is necessary for Vlasov simulation. A nD multiresolution analysis can also be constructed from the 1D multiresolution analysis by tensor products

V j = V j × V j • • • × V j .
Defining Ibin = {I = (i 1 , i 2 , .., i n ) ∈ Z n , I = 0, i k ∈ {0, 1} ∀k}, we obtain the formula

V j = V j-1 + I∈Ibin W I j-1 ,
where

W I j-1 = U i 1 j-1 × • • • × U in j-1 with U i j-1 = V j-1 if i = 0, W j-1 if i = 1.
The decomposition of a function f at the scale j 1 can be expressed as

f (X) = K C j 0 K Φ 0,j 0 K (X) + j 1 -1 j=j 0 I∈Ibin K D I,j K Φ I,j K (X),
where

Φ I,j K (X) = φ i 1 ,j k 1 (x 1 ) • φ i 2 ,j k 2 (x 2 ) . . . φ in,j kn (x n ) with φ i,j k = ϕ j k if i = 0, ψ j k if i = 1.

The adaptive algorithm

In the initialization phase, we first compute the wavelet decomposition of the initial condition f 0 which is known analytically, and then proceed by compressing it, i.e. eliminating the details which are smaller than a threshold that we set. We then construct an adaptive mesh which, from all the possible points at all the levels between our coarsest and finest, contains only those of the coarsest and those corresponding to details which are above the threshold. We denote by G this mesh.

• Prediction in x: We predict the positions of points where the details should be important at the next time split step by advancing in x the characteristics originating from the points of the mesh G. Then we retain the grid points, at one level finer as the starting point, the support of which contains the end point the characteristic.

• Construction of mesh Ĝ: From the predicted mesh G, we construct the mesh Ĝ where the values of the distribution at the next time step shall be computed. This mesh Ĝ contains exactly the points necessary for computing the wavelet transform of f * at the points of G.

• Advection in x: We compute the origin of the characteristics for each point of Ĝ and interpolate its value at the level where it arrives.

• Wavelet transform of f * : We compute the c k and d k coefficients at the points of G from the values of f * at the points of Ĝ.

• Compression: We eliminate the points of G where the details d k are lower than the fixed threshold.

• Computation of the electric field: We compute ρ on the finest grid and solve Poisson.

• Same steps for the v part of the splitting method.

Numerical results

In order to assess the benefits of the adaptive solver we computed the transverse evolution of a semi-Gaussian beam in a uniform and periodic focusing channel. For such a beam the initial distribution function reads

f (r, v) = I πa 2 √ 2πb e -1 2 (v 2 /b 2 ) if r < a,
and f (r, v) = 0 else.

Semi-Gaussian beam in uniform focusing channel

We first consider the evolution of a semi-Gaussian beam in a uniform focusing channel. The tune depression σ/σ 0 is 0,25. Moreover, here, in dimensionless units ω = 2, ω 0 = 1. a = 4/ √ 15, b = 1/(2 * √ 15) et n 0 = 2π. The tests have been performed with a time step of 1/32 th of period, that is ∆t = 2π/(32ω 0 ) = π/16 for a coarse mesh of 8 × 8 points with 5 refinement levels, that is 256 × 256 points on the finest grid, and a refinement threshold of 10 -3 . This means that details such that |d j k | < 10 -3 are considered to be 0 and the corresponding grid point is removed. Figure 2 represents snapshots of the evolution along with the grid points kept by the adaptive algorithm for the computation. We notice that the adaptive grid follows very well the evolution of the fine structures.

Semi-Gaussian beam in periodic focusing channel

We now consider the transverse evolution of a semi-Gaussian beam in a periodic focusing field of the form α(1 + cos 2πz/S) for a tune depression σ/σ 0 of 0.17. Figure 3 represents snapshots of the evolution along with the grid points kept by the adaptive algorithm for the computation. As in the previous case, the adaptive grid follows very well the evolution of the fine structures.

We have introduced in this paper the concept of an adaptive semi-Lagrangian Vlasov solver and proved the feasability of the method for 1D model problems. The method we proposed can be generalized to higher dimensions. Adaptivity is defined through a multiresolution analysis based on interpolating wavelets which is coupled very naturally to the semi-Lagrangian method as fine grid points are added where interpolation on the coarser grid does not do a good enough job.

Adaptive methods enable semi-Lagrangian methods not to waste computing time in regions where nothing is happening and thus make them a lot more efficient. On the other hand, an important overhead is introduced, and programming such a method in an optimal way is a lot more challenging. We are still working on it.

Figure 1 :

 1 Figure 1: Raffinement du maillage en 2D.

Figure 2 :Figure 3 :

 23 Figure 2: Snapshots of the evolution of a semi-Gaussian beam and the associated adatptive grid in a uniform focusing channel

Acknowledgements: This work has been partially supported by the french Commissariat à l'Energie Atomique, Bruyères-Le-Châtel. The authors acknowledge helpful discussions with J.-L. Lemaire.