Sylvain Durand
email: sylvain.durand@mi.parisdescartes.fr

Jalal Fadili
email: jalal.fadili@greyc.ensicaen.fr

Mila Nikolova
email: mila.nikolova@cmla.ens-cachan.fr

Multiplicative noise cleaning via a variational method involving curvelet coefficients

Classical ways to denoise images contaminated with multiplicative noise (e.g. speckle noise) are filtering, statistical (Bayesian) methods, variational methods and methods that convert the multiplicative noise into additive noise (using a logarithmic function) in order to apply a shrinkage estimation for the log-image data and transform back the result using an exponential function. We propose a new method that involves several stages: we apply a reasonable under-optimal hard-thresholding on the curvelet transform of the log-image; the latter is restored using a specialized hybrid variational method combining an ℓ 1 data-fitting to the thresholded coefficients and a Total Variation regularization (TV) in the image domain; the restored image is an exponential of the obtained minimizer, weighted so that the mean of the original image is preserved. The minimization stage is realized using a properly adapted fast Douglas-Rachford splitting. The existence of a minimizer of our specialized criterion and the convergence of the minimization scheme are proved. The obtained numerical results outperform the main alternative methods.

Introduction

In many active imaging systems (e.g. synthetic aperture radar, laser or ultrasound imaging), the data for the unknown image S 0 : Ω → I R + , Ω ⊂ I R 2 , are severely corrupted with multiplicative noise. Then several independent measurements for the same image are needed:

S k = S 0 η k + n k , ∀k ∈ {1, • • • , K}, (1)
where η k : Ω → I R + , and n k represent the multiplicative and a typically zeromean additive noise, ∀k. Commonly (see e.g. [START_REF] Ulaby | Handbook of Radar Scattering Statistics for Terrain[END_REF]) η k is modeled as a onesided exponential probability density function (pdf) (cf. Fig. 1(a)): pdf(η k ) = µ e -µη k 1l I R+ (η k ) for µ = 1. In practice, one takes an average of all measurements, see e.g. Fig. 2(b). Since 1 K K k=1 n k ≈ 0, the data read (cf. e.g. [START_REF] Achim | Sar image denoising via bayesian wavelet shrinkage based on heavy-tailed modeling[END_REF][START_REF] Ulaby | Handbook of Radar Scattering Statistics for Terrain[END_REF][START_REF] Xie | SAR speckle reduction using wavelet denoising and markov random field modeling[END_REF]):

S = 1 K K k=1 S k = S 0 1 K K k=1 η k = S 0 η . (2)
Usually all η k are independent. Denoting by Γ the usual Gamma-function, the mean of the noise η in (2) has a Gamma distribution (cf. Fig. 1(b)):

η = 1 K K k=1 η k : pdf(η) = K K η K-1 exp (-Kη) Γ (K) . (3)
Various adaptive filters have been proposed, see e.g. [START_REF] Krissian | Oriented speckle reducing anisotropic diffusion[END_REF][START_REF] Yu | Speckle reducing anisotropic diffusion[END_REF]: they work well when the noise is moderate or weak, i.e. for K large. Bayesian, variational or diffusionbased methods have been proposed as well; see e.g. [START_REF] Aubert | A variational approach to remove multiplicative noise[END_REF][START_REF] Ma | Combined Curvelet Shrinkage and Nonlinear Anisotropic Diffusion IEEE Trans[END_REF][START_REF] Rudin | Multiplicative denoising and deblurring: Theory and algorithms[END_REF][START_REF] Walessa | Model-based despeckling and information extraction from sar images[END_REF]. Numerous methods convert the multiplicative noise into additive noise by

v = log S = log S 0 + log η = u 0 + n, (4)
see e.g. [START_REF] Achim | Sar image denoising via bayesian wavelet shrinkage based on heavy-tailed modeling[END_REF][START_REF] Fukuda | Suppression of speckle in synthetic aperture radar images using wavelet[END_REF][START_REF] Pizurica | A review of wavelet denoising in mri and ultrasound brain imaging[END_REF][START_REF] Xie | SAR speckle reduction using wavelet denoising and markov random field modeling[END_REF]. Then the pdf of n reads (cf. Fig. 1(c)):

n = log η : pdf(n) = K K Γ (K) -1 exp -K(n -e n ) . (5)
One can prove that E [n] = ψ 0 (K)log K and Var [n] = ψ 1 (K), where ψ k (z) = d dz k+1 log Γ (z) is the polygamma function. A common strategy is to decompose the log-data v into a multiscale frame for L 2 (I R 2 ) (an over-complete basis), say W ≡ {w i : i ∈ I} where I is a set of indexes:

y = W v = W u 0 + W n. (6)
By the Central Limit Theorem, the noise W n in y is nearly Gaussian -cf. Fig.

0 10 1 0 1 2 1 -2 0 1 1 -1 0 1 1 (a) η k (b) η = 1 K P K k=1 η k (c) n = log η (d) W n

1(d).

Then coefficients y are denoised using shrinkage estimators T : I R → I R:

y T [i] = T (W v)[i] , ∀i ∈ I. (7)
Shrinkage functions designed for multiplicative noise were proposed e.g. in [START_REF] Achim | Sar image denoising via bayesian wavelet shrinkage based on heavy-tailed modeling[END_REF][START_REF] Xie | SAR speckle reduction using wavelet denoising and markov random field modeling[END_REF].

Let W ≡ { w i : i ∈ I} be a left inverse of W . Then a denoised log-image v T reads v T = i∈I T ((W v)[i]) w i = i∈I T (y[i]) w i . (8)
Then the sought-after image is of the form S T = exp v T .

Our approach. We apply (4) and consider a tight-frame transform of the logdata. The restored log-image (section 2) minimizes a criterion composed of an ℓ 1 -fitting to the (suboptimally) hard-thresholded frame coefficients and a Total Variation (TV) regularization in the image domain. The minimization (section 3) uses a specialized Douglas-Rachford splitting. The full algorithm, involving a bias correction, is given in section 4. Experiments are presented in section 5. Some notations. (. T ) means transposed, (. * ) means convex conjugate and (. ⋆ ) means adjoint.

Restoration of the log-image

Here we consider how to restore a good log-image given data v : Ω → I R obtained using (4). We focus on methods which, for a given preprocessed data set, lead to convex optimization problems. We comment only variational methods and shrinkage estimators since they underly our specialized hybrid objective function.

Drawbacks of shrinkage restoration and variational methods

Shrinkage restoration. The main problems with these methods, sketched in ( 7)- [START_REF] Chesneau | Stein block thresholding for image denoising[END_REF], is that shrinking large coefficients entails an erosion of the spiky features, while shrinking small coefficients yields Gibbs-like oscillations in the vicinity of edges and a loss of details in the textured area. On the other hand, if shrinkage is insufficient, some coefficients bearing mainly noise can remain almost unchanged-we call such coefficients outliers-and [START_REF] Chesneau | Stein block thresholding for image denoising[END_REF] shows that they yield artifacts with the shape of the functions w i , see Fig. 2. Even though various improvements were brought, these artifacts remain visible-see the results on Fig. Variational methods. In these methods, the restored function minimizes a criterion F v of the form

F v (u) = ρ Ω ψ u(t), v(t) dt + Ω ϕ(|∇u(t)|) dt, (9)
where ψ : I R + → I R + measures closeness to data and ϕ(|∇u(•)|) introduces priors via a trade-off parameter ρ > 0. A classical choice is

ψ = u(•)-v(•) 2 .
It is usually required that the potential function ϕ : I R + → I R + promotes images involving edges. Analysing the minimizers of F v as solutions of PDE's on Ω, Rudin, Osher and Fatemi [START_REF] Rudin | Nonlinear total variation based noise removal algorithm[END_REF] exhibited that ϕ(|∇u(t)|) = |∇u(t)|, leads to such images, where

for any z(t) = (z 1 (t), z 2 (t)) ∈ I R 2 , t ∈ Ω, one sets |z(t)| def = z 1 (t) 2 + z 2 (t) 2 .
The resulting regularization term is known as Total Variation (TV) and will be denoted by • TV . However, whatever smooth data-fitting is chosen, this regularization yields images containing numerous constant regions (called staircasing effect), hence textures and fine details are removed, see [START_REF] Nikolova | Weakly constrained minimization. Application to the estimation of images and signals involving constant regions[END_REF]. The method in [START_REF] Aubert | A variational approach to remove multiplicative noise[END_REF] is of this kind and operates in the image domain; the fitting term is derived from (3) and the denoised image Ŝ, defined by

Ŝ = arg min Σ F S for F S (Σ) = ρ(K) log Σ(t) + S(t)/Σ(t) dt + Σ TV , (10)
exhibits constant regions (see section 5). In [START_REF] Shi | A nonlinear inverse scale space method for a convex mutiplicative noise model[END_REF], the regularization Σ TV is changed into log Σ TV so as to reformulate the model as a convex problem, and not to over smooth the image parts with higher gray values. To recover the denoised image, we applied Ŝ ∝ exp(û) for û = arg min

u where F v (u) = ρ u -v 2 + u TV . (11)
Following [START_REF] Rudin | Nonlinear total variation based noise removal algorithm[END_REF], various edge-preserving convex functions ϕ have been proposed; see [START_REF] Aubert | Mathematical problems in image processing[END_REF] for a recent overview. Even though ϕ ′ (0) = 0 alleviates stair-casing, a systematic drawback of the resulting restored images is that the amplitude of edges is underestimated; thus neat edges or spiky areas are subjected to erosion.

Hybrid methods

Hybrid methods, see e.g. [START_REF] Candès | New multiscale transforms, minimum total variation synthesis. Applications to edge-preserving image reconstruction[END_REF][START_REF] Coifman | Combining the calculus of variations and wavelets for image enhancement[END_REF][START_REF] Durand | Reconstruction of wavelet coefficients using total variation minimization[END_REF][START_REF] Malgouyres | Mathematical analysis of a model which combines total variation and wavelet for image restoration[END_REF], combine the information contained in the large coefficients y[i] obtained according to [START_REF] Chambolle | An algorithm for total variation minimization and application[END_REF] with priors directly on the image u. They amount to define the restored function û by

minimize Φ(u) subject to û ∈ {u : |(W (u -v)) [i]| ≤ µ i , ∀i ∈ I} .
Using an edge-preserving regularization, such as Φ = TV is a pertinent choice. The selection of parameters {µ i } i∈J is more tricky. This choice must take into account the magnitude of the relevant data coefficient y[i]. However, choosing µ i based solely on y[i], as done in these papers, is too rigid since there are either correct data coefficients that incur smoothing (µ i > 0), or noisy coefficients that are left unchanged (µ i = 0). A good compromise that we adopt is to determine (µ i ) i∈I based both on the data and on the prior term.

A specialized hybrid criterion

Given the log-data v obtained by (4), we apply a frame transform as in [START_REF] Chambolle | An algorithm for total variation minimization and application[END_REF] to get y = W v = W u 0 + W n. The noise contained in the i-th datum reads n, w i . The low frequency approximation coefficients carry important information on the image. Therefore, a good choice is to keep them intact at this stage. Let I * ⊂ I denote the subset of all such elements of the frame. Then we apply a hard-thresholding operator T H [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF] to all coefficients I \ I * :

y TH [i] def = T H y[i] , ∀i ∈ I \ I * , where T H (t) = 0 if |t| ≤ T, t otherwise, ( 12
)
where T is an underoptimal threshold in order to preserve the information relevant to edges and to some fine details in textured areas, contained in the small coefficients. Let us consider

v TH = i∈I1 W v[i] w i , where I 1 = {i ∈ I : |y[i]| > T } ∪ I * . ( 13
)
The image v TH contains a lot of artifacts with the shape of the w i for those y[i] that are noisy but above the threshold T , as well as information on the fine details in the original log-image u 0 . In all cases, whatever the choice of T , an image of the form v TH is unsatisfactory-see Fig. 2.

The denoised coefficients, denoted by x, are obtained based on the under-thresholded data y TH . We focus on hybrid methods of the form: x = arg min x F (x) for

F (x) = Ψ (x, y TH ) + Φ( W x)
, where Ψ is a data-fitting term in the frame domain and Φ is an edge-preserving regularization term in the log-image domain. Let us denote

I 0 = I \ I 1 = {i ∈ I \ I * : |y[i]| ≤ T }. (14)
Coefficients y[i] for i ∈ I 0 can be of the two types.

[i] = y[i],
then û would contain an artifact with the shape of w i since by [START_REF] Durand | Denoising of frame coefficients using l1 data-fidelity term and edge-preserving regularization[END_REF] we get v TH = j\i x[j] w j + y[i] w i . Instead, x[i] must be restored according to the prior Φ. This analysis clearly defines the goals that the minimizer x of F is expected to achieve: restored coefficients x[i] have to fit y TH [i] exactly if they are coherent with the prior Φ, otherwise they have to be restored according to Φ. Since [START_REF] Nikolova | Minimizers of cost-functions involving nonsmooth data-fidelity terms. Application to the processing of outliers[END_REF] it is known that such requirements can be satisfied by criteria F where Ψ is non-smooth at the origin (e.g. ℓ 1 ), see also [START_REF] Durand | Denoising of frame coefficients using l1 data-fidelity term and edge-preserving regularization[END_REF]. For these reasons, we focus on

F (x) = Ψ (x) + Φ(x), (15)
where, for Λ = diag(λ i ) i∈I ,

Ψ (x) = i∈I1∪I * λ i |(x -y)[i]| + i∈I0 λ i |x[i]| = Λ(x -y TH ) 1 , (16)
Φ(x) = Ω |∇ W x| ds = W x TV . (17)
In the pre-processing step [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF] we do not recommend the use of a shrinkage function other than T H since it will alter all the data coefficients without restoring them faithfully. Via T H , we base our restoration on data y TH where all nonthresholded coefficients keep the original information on the sought-after image. The theorem stated next addresses the existence and the uniqueness of a minimizer for F . Given y, let G y be the (convex) set of all minimizers of F :

G y def = x ∈ ℓ 2 (I) : F (x) = min x∈ℓ 2 (I) F (x) . ( 18
)
Theorem 1.

[13] For y ∈ ℓ 2 (I) and T > 0 given, consider F as defined in [START_REF] Durand | Multiplicative noise removal using L1 fidelity on frame coefficients[END_REF], where Ω ∈ I R 2 is open, bounded and its boundary ∂Ω is Lipschitz. Suppose that {w i } i∈I is a frame of L 2 (Ω) and the operator W is the pseudo-inverse of W . Assume also that λ min = min i∈I λ i > 0. Then G y is nonempty, and for all

x1 , x2 ∈ G y , ∇ W x1 ∝ ∇ W x2 , a.e. on Ω.
In words, Ŝ1 = W x1 and Ŝ2 = W x2 have the same level lines, i.e. they differ by a local change of contrast; the latter is usually invisible to the naked eye. The choice of λ i is investigated in [START_REF] Durand | Denoising of frame coefficients using l1 data-fidelity term and edge-preserving regularization[END_REF]. Following this analysis, we use only two values for λ i , depending only on the set I ǫ the index i belongs to. We focus on curvelets transforms of the log-data because (a) such a transform captures efficiently the main features of the data and (b) it is a tight-frame which is helpful for the subsequent numerical stage.

Minimization for the log-image

Let Γ 0 (H) denote the class of proper lower-semicontinuous convex functions on a Hilbert space H. Now we focus on the minimization problem find x such that F (x) = min

x

F for F = Ψ + Φ, (19)
where Ψ and Φ are defined in ( 16)-( 17). Clearly, Ψ, Φ ∈ Γ 0 (ℓ 2 (I)), hence F ∈ Γ 0 (ℓ 2 (I)). The set G y in ( 18) is non-empty by Theorem 1 and can be rewritten as G y = {x ∈ ℓ 2 (I) x ∈ (∂F ) -1 (0)}, where ∂F stands for subdifferential.

Minimizing F amounts to finding a solution to the fixed point equation

x = (Id + γ∂F ) -1 (x) , (20)
where (Id + γ∂F ) -1 is the resolvent operator associated to ∂F , γ > 0 is the proximal stepsize and Id is the identity map on ℓ 2 (I). Since (Id + γ(∂Ψ + ∂Φ)) -1 cannot be calculated in closed-form, we focus on splitting methods that use separately the resolvent operators (Id + γ∂Ψ ) -1 and (Id + γ∂Φ)) -1 .

Specialized Douglas-Rachford (D-R) splitting algorithm

The D-R family is the most general class of monotone operator splitting methods. Given a sequence µ t ∈ (0, 2), D-R methods can be expressed via the recursion

x (t+1) = 1- µ t 2 Id+ µ t 2 (2(Id+γ∂Ψ ) -1 -Id) • (2(Id+γ∂Φ) -1 -Id) x (t) . ( 21
)
Since problem [START_REF] Malgouyres | Mathematical analysis of a model which combines total variation and wavelet for image restoration[END_REF] has solutions, we have the following convergence result:

Theorem 2. Let γ > 0 and µ t ∈ (0, 2) be such that t∈I N µ t (2 -µ t ) = +∞. Take x (0) ∈ ℓ 2 (I) and consider the sequence of iterates defined by [START_REF] Nikolova | Minimizers of cost-functions involving nonsmooth data-fidelity terms. Application to the processing of outliers[END_REF]. Then, (x (t) ) t∈I N converges weakly to some point x ∈ ℓ 2 (I) and (Id+γ∂Φ) -1 (x) ∈ G y .

The statement follows from [10, Corollary 5.2]. The sequence µ t = 1, ∀t ∈ I N fits.

Proximal calculus

Proximity operators, invented in [START_REF] Moreau | Fonctions convexes duales et points proximaux dans un espace hilbertien[END_REF], generalize convex projection.

Definition 1 (Moreau [20]). Let ϕ ∈ Γ 0 (H). Then ∀x ∈ H the function z → ϕ(z)+ x-z 2 /2, for z ∈ H, achieves its infimum at a unique point denoted by prox ϕ x. The relevant operator prox ϕ : H → H is the proximity operator of ϕ.

By the minimality condition for prox ϕ , it is easy to see that ∀x, p ∈ H we have p = prox ϕ x ⇐⇒ x -p ∈ ∂ϕ(p) ⇐⇒ (Id + ∂ϕ) -1 = prox ϕ . By introducing the reflection operator rprox ϕ def = 2prox ϕ -Id, the D-R iteration (21) reads

x (t+1) = 1 - µ t 2 Id + µ t 2 rprox γΨ • rprox γΦ x (t) . (22)
Proximity operator of

Ψ Lemma 1. Let x ∈ ℓ 2 (I). Then prox γΨ (x) = y TH [i] + T S γλi (x[i] -y TH [i]) i∈I ,
where

T S γλi (z[i]) = max 0, z[i] -γλ i sign(z[i]) .
The proof is quite standard and can be found in our Report [START_REF] Durand | Multiplicative noise removal using L1 fidelity on frame coefficients[END_REF]. Note that

rprox γΨ (x) = 2 y TH [i] + T S γλi (x[i] -y TH [i]) i∈I -x . (23)
Proximity operator of Φ. Clearly, Φ(x) = • TV • W (x). Computing prox γΦ for an arbitrary W may be intractable. We assume that

(w1) W : ℓ 2 (I) → L 2 (Ω) is surjective; (w2) W W = Id and W = c -1 W ⋆ for 0 < c < ∞; note that W ⋆ W = c Id; (w3) W is bounded. Let X = L 2 (Ω) × L 2 (Ω), •, • X be the inner product in X and • p , p ∈ [1, ∞] the L p -norm on X . Define B γ ∞ (X ) as the γ-radius closed L ∞ -ball in X , B γ ∞ def = z ∈ X : z ∞ ≤ γ = z = (z 1 , z 2 ) ∈ X : |z(t)| ≤ γ, ∀t ∈ Ω , and P B γ ∞ (X ) : X → B γ ∞ (
X ) the associated projector.

Lemma 2. Let x ∈ ℓ 2 (I) and B γ ∞ (X ) is as defined above. Then:

prox γΦ (x) = Id -W • Id -prox c -1 γ • TV • W (x) ; (24)
prox c -1 γ • TV (u) = u -P C (u) , (25)
where

C = div(z) ∈ L 2 (Ω) z ∈ C ∞ c (Ω × Ω), z ∈ B γ/c ∞ (X ) . (26
(z) = 0 if z ∈ C, +∞ if z ∈ C.
where C is given in [START_REF] Shi | A nonlinear inverse scale space method for a convex mutiplicative noise model[END_REF]. Using Definition 1, prox c -1 γ . TV * = P C . Identifying c -1 γ . TV with ϕ and c -1 γ . TV * with ϕ * , equation ( 27) leads to (ii) 4 . ⋄ ¿From ( 24)- [START_REF] Rudin | Nonlinear total variation based noise removal algorithm[END_REF] we easily find that

rprox γΦ (x) = Id -2W • P C • W (x) . ( 28
)
Calculation of the projection P C in (25) on a discrete grid. In this case, W is an M×N tight frame with M= #I ≫ N = #Ω and assumption (w2) reads W W = Id and

W = c -1 W T , c ∈ (0, ∞) hence W T W = c Id). The discrete counterpart of X is X = ℓ 2 (Ω) × ℓ 2 (Ω).
We denote the discrete gradient by ∇ (cf. [START_REF] Chambolle | An algorithm for total variation minimization and application[END_REF] or [START_REF] Welk | Locally analytic schemes: A link between diffusion filtering and wavelets shrinkage[END_REF]) and the discrete divergence Div : X → ℓ 2 (Ω) is defined as Div = -∇⋆ . Moreover, C in (26) admits a simpler expression:

C = Div(z) ∈ ℓ 2 (Ω) z ∈ B γ/c ∞ (X ) , (29)
where

B γ/c
∞ (X ) is defined using the new discrete notations. The projection P C in [START_REF] Rudin | Nonlinear total variation based noise removal algorithm[END_REF] does not admit an explicit form so we provide an iterative scheme for its calculation in the next lemma. Lemma 3. We adapt all assumptions of Lemma 2 to the new discrete setting, as explained above. Consider the forward-backward iteration

z (t+1) = P B 1 ∞ (X ) z (t) + β t ∇ Div(z (t) ) -cu/γ (30) for 0 < inf t β t ≤ sup t β t < 1/4 (31)
where

∀(i, j) ∈ Ω , P B 1 ∞ (X ) (z)[i, j] = z[i, j] if |z[i, j]| ≤ 1; z[i, j]/|z[i, j]| otherwise . ( 32
)
Then (i) (z (t) ) t∈I N converges to a point ẑ ∈ B 1 ∞ (X ); (ii) c -1 γDiv(z (t) ) t∈I N converges to c -1 γDiv(ẑ) = (Id -prox c -1 γ • TV )(u).
The proof of this lemma can be found in our Report [START_REF] Durand | Multiplicative noise removal using L1 fidelity on frame coefficients[END_REF]. The iteration proposed in [START_REF] Xie | SAR speckle reduction using wavelet denoising and markov random field modeling[END_REF] to compute the proximity operator of the TVnorm is different from the projection algorithm of [START_REF] Chambolle | An algorithm for total variation minimization and application[END_REF]. A similar iteration was proposed in [START_REF] Chambolle | Total variation minimization and a class of binary MRF models[END_REF] and in some other articles. The proof we gave is however simpler as it uses known properties of proximity operators. Note that computing prox • TV amounts to solving a discrete ROF-denoising. Our iteration to solve this problem is one possibility among others, see e.g. a recent report [START_REF] Aujol | Some algorithms for total variation based image restoration[END_REF]. A crucial property of the D-R scheme ( 22) is its robustness to numerical errors that may occur when computing the proximity operators prox Ψ and prox Φ , see [START_REF] Combettes | Solving monotone inclusions via compositions of nonexpansive averaged operators[END_REF]. More precisely, let a t ∈ ℓ 2 (I) be an error term that models the inexact computation of prox γΦ in [START_REF] Rudin | Multiplicative denoising and deblurring: Theory and algorithms[END_REF], as the latter is obtained through [START_REF] Xie | SAR speckle reduction using wavelet denoising and markov random field modeling[END_REF]. If the sequence of error terms (a t ) t∈I N and stepsizes (µ t ) t∈I N in Theorem 2 obey t∈I N µ t a t < +∞, then the D-R algorithm ( 22) converges [10, Corollary 6.2]. In our experiments, using 200 inner iterations in [START_REF] Xie | SAR speckle reduction using wavelet denoising and markov random field modeling[END_REF] is sufficient to satisfy this requirement.

Bias correction to recover the sought-after image

Recall from (4) that u 0 = log S 0 and set û = W x(NDR) as the estimator of u 0 , where N DR is the number of D-R iterations in [START_REF] Nikolova | Weakly constrained minimization. Application to the estimation of images and signals involving constant regions[END_REF]. Unfortunately, the estimator û is prone to bias, i.e. E [û] = u 0 -b û. A problem that classically arises in statistical estimation is how to correct such a bias. More importantly is how this bias affects the estimate after applying the inverse transformation, here the exponential. Our goal is then to ensure that for the estimate Ŝ of the image, we have E Ŝ = S 0 . Expanding Ŝ in the neighborhood of E [û], we have

Fig. 1 .

 1 Fig. 1. Noise distributions.

  3(d) and Fig.4(c) in Section 5 using the very recent Stein-block thresholding[START_REF] Chesneau | Stein block thresholding for image denoising[END_REF].

  (a) Noisy, K = 10 (b) T = 2 p Var [n] (c) T = 4 p Var [n] (d) T = 6 p Var [n]

Fig. 2 .

 2 Fig. 2. (a) Noisy Lena obtained according to (1)-(2) for K = 10. (b)-(d) Restorations exp vT H where data v are denoised by hard-thresholding of its curvelet coefficients, see (12)-(13), for different choices of T .

)

  Sketch of the proof. By (w1), range( W ) = L 2 (Ω). Using that domain( • TV ) = L 2 (Ω), we find cone dom • TVrange W = {0}. Statement (i) follows from applying [11, Proposition 11] whose requirements are satisfied. If ϕ ∈ Γ 0 (L 2 (Ω)) and ϕ * is its convex conjugate, the Moreau decomposition [20, Proposition 4.a] asserts prox ϕ + prox ϕ * = Id . (27) Since the conjugate function of a norm is the indicator function of the ball of its dual norm, c -1 γ • TV *

E

  e û = exp (E [û])(1+Var [û] /2+R 2 ) = S 0 exp (-b û)(1+Var [û] /2+R 2 ) , (33)where R 2 is expectation of the Lagrange remainder in the Taylor series. One can observe that the posterior distribution of û is nearly symmetric, hence R 2 ≈ 0. Then b û ≈ log(1v +Var [û] /2) ensures unbiasedness. Consequently, finite sample (nearly) unbiased estimates of u 0 and S 0 are respectively û + log(1 + Var [û] /2), and exp (û) (1 + Var [û] /2). Var [û] can be reasonably estimated by ψ 1 (K), the variance of the noise n in (4) being given in[START_REF] Achim | Sar image denoising via bayesian wavelet shrinkage based on heavy-tailed modeling[END_REF]. Thus, given the restored logimage û, our denoised image read:Ŝ = exp (û) (1 + ψ 1 (K)/2) .(34)(a) Original (256 × 256) (b) Noisy: µ = 1, K = 10 (c)û by (11) and Ŝ by (34) see (2)-(3) psnr=26.2 db, mae=8.5 (d) Stein-block [8] (e) AA algorithm [2] (f) Our method psnr=25.5 db, mae=9.4 psnr=25.4 db, mae=9.4 psnr=26.05 db, mae=8.8

Fig. 3 .

 3 Fig.3. Restoration of (b) using modern methods. Note that (c) is a slightly improved version of[START_REF] Shi | A nonlinear inverse scale space method for a convex mutiplicative noise model[END_REF] and that the restoration in (d) is done in the curvelet domain.

  (a) Original (b) Noisy (c) Stein-block th. (d) Our method

Fig. 4 .

 4 Fig. 4. (a) Shepp-Logan phantom (256 × 256). (b) Noisy, K = 10. (c) Denoised with Stein-block thresholding in the curvelet domain [8] PSNR=24.73dB, MAE=4. (d) Denoised with our algorithm PSNR=31.25dB, MAE=1.87.

  1. Coefficients y[i] bearing mainly noise-then the best choice is x[i] = 0; 2. Coefficients y[i] relevant to edges and other details in u 0 . Since y[i] is difficult to distinguish from the noise, the relevant x[i] should be restored using the edge-preserving prior Φ. Note that a careful restoration must find a nonzero x[i] in order to avoid Gibbs-like oscillations in û. Coefficients y[i] for i ∈ I 1 are of the following two types. 1. Large coefficients which carry the main features of the sought-after function. They verify y[i] ≈ w i , u 0 and can be kept intact. 2. Coefficients highly contaminated by noise, i.e. |y[i]| ≫ | w i , u 0 |. We call them outliers because if we had x

Note that our argument[START_REF] Ulaby | Handbook of Radar Scattering Statistics for Terrain[END_REF] to compute prox c -1 γ • TV (u) is not used in[START_REF] Chambolle | An algorithm for total variation minimization and application[END_REF], which instead uses conjugates and bi-conjugates of the objective function.

Full algorithm to suppress multiplicative noise

Piecing together Lemmas 1 and 2, and Theorem 2, we write down the full multiplicative noise removal algorithm: Task: Denoise an image S corrupted with multiplicative noise according to [START_REF] Aubert | A variational approach to remove multiplicative noise[END_REF]. Parameters: The observed noisy image S, number of iterations N DR (Douglas-Rachford outer iterations) and N FB (Forward-Backward inner iterations), stepsizes µ t ∈ (0, 2), 0 < β t < 1/4 and γ > 0, tight-frame transform W and initial threshold T (e.g. T = 2 ψ 1 (K)), regularization parameters λ 0,1 associated to the sets I 0,1 . Specific operators:

(c) ∇ and Div-the discrete versions of the continuous operators ∇ and div. (d) ψ 1 (•) defined according to (1) (built-in Matlab function).

Output:

Experiments

In all experiments, our algorithm was run using second-generation curvelet tight frame along with the following set of parameters: ∀t, µ t ≡ 1, β t = 0.24, γ = 10 and N DR = 50. The initial threshold T was set to 2 ψ 1 (K). For comparison purposes, some very recent multiplicative noise removal algorithms from the literature are considered: the AA algorithm [START_REF] Aubert | A variational approach to remove multiplicative noise[END_REF] minimizing the criterion in [START_REF] Combettes | Solving monotone inclusions via compositions of nonexpansive averaged operators[END_REF], and the Stein-block denoising method [START_REF] Chesneau | Stein block thresholding for image denoising[END_REF] in the curvelet domain, applied on the log transformed image. The latter is a sophisticated shrinkage-based denoiser that thresholds the coefficients by blocks rather than individually, and has been shown to be nearly minimax over a large class of images in presence of various additive bounded noises. We also tried the L2-TV method where the restored log-image û minimizes [START_REF] Combettes | A Douglas-Rachford splittting approach to nonsmooth convex variational signal recovery[END_REF] and the denoised image Ŝ involves the bias correction (34). Thanks to the bias correction, it can be seen as an improved version of the first method proposed in the recent Report [26, § 4.1]. For fair comparison, the hyperparameters for all competitors were tweaked to reach their best level of performance on each noisy realization.

The denoising algorithms were tested on two images: Lena and Boat, all of size 256×256 and gray-scale in the range [START_REF] Achim | Sar image denoising via bayesian wavelet shrinkage based on heavy-tailed modeling[END_REF]256]. For each image, a noisy observation is generated by multiplying the original image by a realization of noise according to The results are depicted in Figs. 3 and4. Note that the AA algorithm tends to over-regularize the solution. Our denoiser clearly outperforms its competitors.


domain: chim sde

text:
Soulier Coralie

published or not. The documents may come

Analytical and processing methodologies for non-target

screening by coupling passive sampling and high resolution mass spectrometry

Analytical and processing methodologies for non-target screening by coupling passive sampling and high resolution mass spectrometry Soulier Coralie and Togola Anne BRGM Laboratory Division, 3 avenue C. Guillemin, 45100 Orléans, FRANCE E-mail contact: c.soulier@brgm.fr

Nowadays with technological advances, the use of environmental forensic approaches could help to characterize the various sources of groundwater contamination. This implies the need of specific analytical methodology to identify micropollutants, emerging substances or transformation products present at low concentrations. The high resolution mass spectrometry (HRMS) has gained increasingly in importance for monitoring organic compounds. Its high resolving power, mass accuracy and the sensitive full spectrum acquisition are the key points. Contamination profile and pattern of a specific site could be highlighted by this technique with the use of automatic data processing softwares.

The aim is to support public policy development by highlighting and identifying relevant compounds to be monitored in groundwater. The main difficulties for the implementation of monitoring are sometimes low and fluctuating concentration levels and complex mixture of pollutants. No therefore there is a strong interest to combine passive sampling to HRMS. Passive samplers allow accumulating compounds during exposure that improve trace detection and integrating pollution fluctuations. The Polar Organic Chemical Integrative Sampler (POCIS) was employed to sampling polar and semi-polar compounds (pesticides, pharmaceuticals, phenolic compounds, triazoles….).

Different sites impacted by agricultural, urban or industrial pollution sources were investigated and sampled during several months. Grab and passive sampling were deployed and analyzed by LC-QTof. To process data, different approaches were investigated according to identification confidence of NORMAN Network (Schymanski et al., 2014). The first one is based on research from target compounds listed on our homemade database (around 450 with experimental data on our system as retention time, exact masses for molecular and fragment ions) or in bibliography. The second one concerns non-targeted screening. The identification of unknown compounds was tried by using several databases such as Norman Mass Bank or Chemspider and different application manager for MassLynx software.

This approach allows highlighting the use of passive samplers as storage tool becausemore compounds are identified with POCIS. Transformation products of pesticides or pharmaceutical compounds were present in most of samples.

Keywords: non-target screening, passive sampling, POCIS, high resolution mass spectrometry Bibliography:

Schymanski E., Jeon J., Gulde R., Fenner K., Ruff M., Singer H. and Hollender J.( 2014), Identifying small molecules via high resolution mass spectrometry: communicating confidence. Environmental Science andTechnology, 48, 2097-2098.


H Aochi

T Ulrich

D Michea

F Dupros

S06b Strong Ground Motion: Earthquake Scenarios

Keywords:

Variability of ground motion estimation according to various earthquake scenarios is studied in the region of the Sea of Marmara around Istanbul, Turkey, through numerical simulations using a 3D finite-difference method with CPU-GPU implementation. The earthquake scenarios dynamically simulated based on mechanical models are given with probability according to the geological, geophysical and seismological hypotheses. We adopt a 3D geological model of the Sea of Marmara taken into account of the basin structure (down to about 10 km depth) and the existence of the Sea layer (1.5 km depth at max.). Although beneath the land part around the Sea of Marmara is briefly a 1D structure, the influence of the 3D structure of the basin on the ground motion caused by the earthquakes along the North Anatolian fault is significant. Various earthquake scenarios of a magnitude of about 7 lead to a very heterogeneous ground motion pattern. For the selected points, we carry out statistical analysis on its characterization. In Istanbul area, a very strong ground motion (peak ground velocity larger than 1 m/s) may not be the most probable scenario, but has a probability of a few percents in the case where the stress drop is extremely high enough to lead to a super-shear rupture propagation.


EA Philippe Martin Clillac-Arp
email: philippe.martin@linguist.jussieu.fr

The Autosegmental-Metrical Prosodic Structure: not fit for French?

Keywords: sentence intonation, autosegmental-metrical, ToBI, prosodic structure, French

Since two decades at least, intonation phonology is almost exclusively dominated by the Autosegmental-Metrical model, which appears universally appropriate to describe sentence intonation in most languages. However, careful examination of the AM limitations and drawbacks may question this universality, leading way to improve and possibly modify key properties of this theoretical approach, at least for languages such as French.

Introduction

Since two decades at least, the so-called Autosegmental-Metrical (AM) model [START_REF] Selkirk | On prosodic structure and its relation to syntactic structure[END_REF] has been dominant in intonation phonology. In this model, the prosodic structure organizes hierarchically prosodic events (PE) in three non-recursive levels: a first level assembles syllables σ, content words Wc (verbs, nouns adjectives and adverbs) and function words Wf (conjunctions, pronouns,…) into accentual phrases (AP); a second level groups AP into intonation phrases (IP) (Fig. 1); finally a phonological utterance (PU) eventually groups sequences of IP.

Figure 1. Autosegmental-Metrical prosodic structure [START_REF] Selkirk | On prosodic structure and its relation to syntactic structure[END_REF] The prosodic events PE are aligned on accentual phrases specific syllables and are described as sequences of tones belonging to the ToBI notational system (tones and break indices). This system uses High (H) and Low (L) symbols to transcribe melodic targets as perceived or observed on fundamental frequency curves obtained from the speech signal acoustic analysis.

Since the AM approach is widely used on many languages other than English (Japanese [START_REF] Pierrehumbert | Japanese Tone Structure[END_REF], Korean [START_REF] Jun | The Phonetics and Phonology of Korean Prosody: intonational phonology and prosodic strcture[END_REF], Dutch [START_REF] Gussenhoven | Transcription of Dutch Intonation[END_REF], German [START_REF] Baumann | GToBI -a phonological system for the transcription of German intonation[END_REF], Italian [START_REF] Avesani | ToBIt: un sistema di trascrizione per l'intonazione italiana[END_REF], French [START_REF] Sun-Ah | The Realizations of the Accentual Phrase in French Intonation[END_REF], etc.), one may think that it carries some universal value by capturing in its essence basic characteristics of sentence intonation of so many languages. However, open or allusive critics have been emerging for some time (e.g. [START_REF] Ladd | Intonational phrasing: The case for recursive prosodic structure[END_REF]) leading to a more serious analysis of the AM limitations and drawbacks, and these critics may apply particularly to French.

Anything wrong with the AM Prosodic Structure?

Indeed, since the AM approach is so widely used in intonation sentence phonology and applied to many different languages, what could be wrong with it? I will try to enumerate a few possibly questionable points:

a. The AM prosodic structure is non-recursive. This property has already been discussed by various authors, and among other reasons originates in my opinion from the fact that very short sentences were used as experimental justification for this property. Indeed, exclusive use of limited length sentences prevents the observation of recursiveness, even in English [START_REF] Ladd | Intonational Phonology[END_REF];

b. Descriptions of prosodic events underlying a prosodic structure do not take duration parameters into account. The ToBI system has no explicit provision to describe temporal aspects of sentence intonation other that the perceived break durations (which is seldom used); c. While other transcription systems are either available or could be more or less easily adapted to fit specific properties of a given language, the quasi exclusive use of the ToBI system involves an oversimplification of the description of melodic events, oversimplification sometimes compensated at a later stage by complex tone alignment rules aimed to better take the phonetic details of melodic movements into account; d. In many instances, confusion exists between phonological and phonetic descriptions of prosodic events. Some authors give PE descriptions so detailed that the appear purely phonetic rather than phonological (e.g. [START_REF] D'imperio | Italian intonation: An overview and some questions[END_REF]); e. Contextual properties of prosodic events are often ignored, as there seem to be a strong underlying assumption that prosodic events share properties similar to phonemes. This aspect is intriguing, as the AM approach was proposed to address the possible effect of context in the realization of melodic contours. Contextual rules may appear only on the surface structure, but are generally not considered as an inherent property of PE;

f. In early versions of the AM framework, the prosodic structure was assumed to be congruent with the sentence syntactic structure. This implies that only one prosodic structure could be associated to a given sentence. Even if congruence with syntax is not necessarily retained today as an hypothesis, it is rare to find an author considering the possibility to associate more than one prosodic structure to a given syntactic structure. A notable exception can be found in [START_REF] Post | Restructured phonological phrases in French: Evidence from clash resolution[END_REF]; g. As other less known theoretical approaches, AM ignores a basic property of sentence intonation, i.e. to be encoded by prosodic events encoded and decoded sequentially by the speaker and the listener. Therefore, it may be misleading to consider prosodic events on a piece of paper as emerging at once to represent the prosodic structure, as they appear in reality in a timely fashion one after the other. This time domain dynamic aspect may modify the way we envision sentence intonation and the prosodic structure [START_REF] Martin | Intonation du français[END_REF];

h. A last point pertains to the quasi exclusive use of laboratory speech, generally involving (very) short sentences. This limited choice of data, justified in the early years by technical limitations, prevents the AM prosodists to observe data that would seriously question the use of their approach.

I will now discuss in some details these points, while trying to suggest possible modifications and improvements.

AM intonation phonology and French experimental data

The AM prosodic structure is non-recursive

Critics of the non-recursiveness constrain appear as early as 1986 (see Ladd [START_REF] Ladd | Intonational phrasing: The case for recursive prosodic structure[END_REF]), but this author later revised his opinion on the recursiveness of the prosodic structure [START_REF] Ladd | Intonational Phonology[END_REF]. In a recent work on French sentence intonation [START_REF] Michelas | Durational cues and prosodic phrasing in French: evidence for the intermediate phrase[END_REF], despite the use of very short sentences in the experiment presented in this study, the need for at least an intermediate IP appears clearly. This point can be also shown by comparing the two following examples in French:

a) Le marin roumain n'aurait pas voulu ranimer la jolie maman.

(The Romanian sailor would not have wanted to revive the beautiful mother).

b) Si le rat marron avait voulu manger le long mulot le marin roumain n'aurait pas voulu ranimer la jolie maman.

(If the brown rat had wanted to eat the long field mouse the Romanian sailor would not have wanted to revive the beautiful mother).

By adding the relative phrase Si le rat marron avait voulu manger le long mulot in front of sentence a) reveals the need to consider at least one extra level in the prosodic structure, Indeed, si le rat marron and avait voulu manger le long mulot appears as second level IP's grouped into a first level IP Si le rat marron avait voulu manger le long mulot. Prosodically, marron bears a falling melodic contour whereas the contour on roumain is rising, indicating the existence of two IP's (Fig. 4).

The ToBi transcription does not transcribe duration

Since ToBI notation system uses only combinations of High and Low tones, no provision is made to represent possible contrasts in duration or in amplitude of melodic variation that would differentiate prosodic events. Rhythmic factors are also absent from the transcription, except indirectly through the scale of breaks (which are based on perception). Fig. 5 shows how two melodic contours of different durations and melodic variation can be transcribed by the same sequence of tones.

ToBI use with the AM prosodic structure tend to confuse phonetic and phonology

Since the AM prosodic structure has by definition only two levels, only three phonologically distinct prosodic events should a priori be considered to encode this structure. However, in many studies, phonetic details of the melodic curve are given by defining the position of the High or Low target tone inside the stressed syllable.

The transcription of the prosodic event on the stressed syllable of marron and roumain in the example of Fig. 6 can be phonological by using H* in both cases but aligned on the beginning of the syllable for marron, and H* aligned on the end of the syllable for roumain. In either case, the alignment of the target H or L can give a proper account to either a rising or falling melodic contour by selecting an adequate phonetic alignment on the left or the right syllable boundary. Although this can describe data properly, it fails to capture the contrast of melodic slope, characteristic of French intonation [START_REF] Martin | Intonation du français[END_REF]. (eurhythmic alignment). More than one prosodic phrasing are then possible for this example (Fig. 7).

Figure 7. Two prosodic structures that can be associated with the same sentence "Marion adore le whisky écossais".

The AM approach does not provide for the dynamic time aspects of the prosodic structure

In the AM framework, the sentence prosodic structure is viewed globally, taking into account all the prosodic events at once from the beginning to the end, without taking into account the sequence of events in function of time. However, from the point of view of the speaker and the listener, the situation is quite different: whereas the speaker can achieve some planning ahead in the production of the prosodic structure, it is barely the case for the listener, who has to process the linguistic information from the sequence of units perceived one by one along the time scale. In this process, prosodic events are used as signals triggering partial processing of the already perceived syllables, by concatenation of strings of already stored units to form larger syllabic groups organized in stress groups. This process is totally absent in the AM framework.

Examples of laboratory speech given in the AM literature are too short

The apparent justification to limit to two the number of levels of the AM prosodic structure may stem from the fact that laboratory speech analysis of sentence intonation was performed on short sentences with a limited syntactic complexity.

In French for instance, examples analyzed in AM studies could be as simple as Marion mangera des bananes ("Marion will eat bananas") [START_REF] Jun | The Phonetics and Phonology of Korean Prosody: intonational phonology and prosodic strcture[END_REF], Le mari d'Amanda réclamait sa bicyclette ("Amanda's husband reclaimed his bicycle") or Marie a rencontré les amis de Rémy dans la journée ("Mary met Rémy's friends in the afternoon") [START_REF] Michelas | Durational cues and prosodic phrasing in French: evidence for the intermediate phrase[END_REF]. Likewise, in Italian, complex studies were carried on sentences like Mamma andava a ballare da Lalla or Io dicevo mamma [START_REF] D'imperio | Italian intonation: An overview and some questions[END_REF]. In Spanish, Cuando hubo hablado, Juan se fue or Le dieron el número de vuelo [START_REF] Mary | Intonation across Spanish, in the Tones and Break Indices framework[END_REF].

The analysis of this kind of examples obviously prevents the observation of a prosodic structure with more than one level IP. It is time for the AM approach to confront the theory to more realistic data.

Transcription as theory

Badiou [START_REF] Badiou | Le concept de modèle -Introduction à une épistémologie matérialiste des mathématiques[END_REF] and Ochs [START_REF] Ochs | Transcription as theory[END_REF] demonstrated that the choice of a transcription system determines the theory that uses this system, whereas, in linguistics in particular, it should be the reverse: the models derived from a theory should determine the transcription system to analyze the experimental data.

It appears rather obviously that the ToBI transcription system acts as a filter retaining only selected parts of the information present in the data. For example when for a language such as French, prosodic events may contrast by the range of variation of melodic contour or by a concave shape; these differences cannot be easily captured with a combination of the available symbols H and L and their variants.

By specifying the alignment of the tonal target inside the syllable, it is possible, although in a not very intuitive way, to give an account for complex Fo curved shapes. A recent encounter to define a French ToBI set of contours held in Tarragona [START_REF]French ToBI, Workshop on Romance ToBI, Satellite workshop PaPI[END_REF] revealed the difficulty to obtain satisfactory set of prosodic events description for French using ToBI symbols.

To summarize, the use of ToBI system not only filters the data, but obscures the explanation principles underlying the concept of prosodic structure applied to the sentence.

The AM approach is short in explanation principle

The AM prosodic structure lacks a convincing explanation principle: what it tells us phonologically is that intonation sentence assembles prosodic words (equivalent to AP) in two levels. Therefore, the only possible falsification test to ensure the AM approach is not tautological consists to demonstrate that indeed only two levels are found in the prosodic structure of any language.

Furthermore, no explanation is proposed pertaining to the phonological role of the melodic contours as described by ToBI notation in the indication of the structure. Prosodic words and their contours actually appear independent from the structure that organize their hierarchy in the sentence, whereas it may seem intuitively (and demonstrated by careful data analysis) that indeed melodic contour act as indicator of the prosodic structure [START_REF] Martin | Intonation du français[END_REF]. Furthermore, when the prosodic structure becomes more complex, the contrasts between the realizations of melodic contours become more elaborate and involve more phonetic features, an easily observable fact that is not taken into account by the AM approach.

Conclusions

The prosodic structure as defined in the Autosegmentalmetrical approach operates as a filter allowing the phonologist and possibly the phonetician to select pertinent characteristics from the complexity of prosodic data, fundamental frequency, duration and intensity as revealed by acoustical analysis. As such, it allowed a new generation of researchers to enter the relatively new field of prosody offered by the development of computer technology and the availability of acoustic analysis software such as Praat or WinPitch.

Still the time has come to go beyond the AM approach to take into account new insights given by a critical analysis of the properties and limitations of the AM prosodic structure. In particular, the point of view of the listener and the relative independence of syntax as revealed by spontaneous speech data should definitely enter in a new theoretical approach of sentence prosody and phrasing.

Figure 2 .

 2 Figure 2. Melodic curve of the French example « Si le rat marron avait voulu manger le long mulot le marin roumain n'aurait pas voulu ranimer la jolie maman » with the stressed syllables between vertical cursors highlighted.

Figure 3 .

 3 Figure3. The AM prosodic structure of the example of Fig.2with both IP's at the same level.

Figure 4 .

 4 Figure 4. The prosodic structure of the example of Fig. 2 as indicated by the melodic contours system of contrasts. Extra levels are needed to adequately represent the structure encoded by the melodic contours.

Figure 5 .

 5 Figure 5. The same L H* ToBI tone sequences transcribe contours with different durations.

Figure 6 . 3 . 4 .

 634 Figure 6. The same L H* ToBI tone sequences transcribe contours with rising or falling melodic slope. 3.4. In the AM approach only one prosodic structure can be associated with a given syntactic structure In most if not all AM studies it is assumed that only one prosodic structure can be assigned to a given text (or a given syntactic structure). Still factors like eurhythmicity are clearly interacting in the speaker realization, and this is particularly obvious in French. In Marion adore le whisky écossais ("Marion loves Scotch whisky") for example, prosodic phrasing can be [Marion] [adore le whisky écossais] (syntactic alignment) as well as [Marion adore] [le whisky écossais]


Pratap Tokekar

Volkan Isler

Antonio Franchi

Multi-Target Visual Tracking With Aerial Robots

We study the problem of tracking mobile targets using a team of aerial robots. Each robot carries a camera to detect targets moving on the ground. The overall goal is to plan for the trajectories of the robots in order to track the most number of targets, and accurately estimate the target locations using the images. The two objectives can conflict since a robot may fly to a higher altitude and potentially cover a larger number of targets at the expense of accuracy.

We start by showing that k ≥ 3 robots may not be able to track all n targets while maintaining a constant factor approximation of the optimal quality of tracking at all times. Next, we study the problem of choosing robot trajectories to maximize either the number of targets tracked or the quality of tracking. We formulate this problem as the weighted version of a combinatorial optimization problem known as the Maximum Group Coverage (MGC) problem. A greedy algorithm yields a 1/2 approximation for the weighted MGC problem. Finally, we evaluate the algorithm and the sensing model through simulations and preliminary experiments.

I. INTRODUCTION

We study the problem of tracking multiple moving targets using aerial robots. We consider the scenario where cameras that face downwards are mounted on the robots to track targets moving on the ground plane. A robot can potentially track more targets by flying to a higher altitude, thus increasing its camera footprint. However, this may reduce the quality of the view due to the increased distance between the cameras and the targets. There is a trade-off between the number of targets tracked and the corresponding quality of tracking. We investigate this trade-off and present an approximation algorithm for multi-target tracking.

We start by showing that it may not always be possible to track all targets while always maintaining the optimal quality of tracking (or any factor of the optimal quality), even if the targets' motion is fully known. Hence, we focus on the following two variants: maximize the number of targets tracked subject to a desired tracking quality per target, and maximize the sum of quality of tracking for all targets. The two problems can be formulated as the unweighted and weighted versions of the Maximum Group Coverage Problem (MGC). A simple greedy approach provides a 1/2 approximation to unweighted MGC [START_REF] Chekuri | Maximum coverage problem with group budget constraints and applications[END_REF]. We show that the approximation guarantee also holds for the weighted case which allows a practical solution to the trajectory planning problem with provable performance guarantees. We evaluate the algorithm in simulations and preliminary experiments with an indoor platform using four aerial robots.

The rest of the paper is organized as follows. We begin with the related work in Section II. The problem setup and a discussion of the sensing quality are presented in Section III. The infeasibility of tracking all targets with a constant factor of the optimal quality is proven in Section IV. The tracking algorithm is presented in Section V, and evaluated through simulations and preliminary experiments in Sections VI and VII respectively. Section VIII concludes the paper.

II. RELATED WORK

Target tracking is an important problem for robotics, and has been widely studied under different settings. Spletzer and Taylor [START_REF] Spletzer | Dynamic sensor planning and control for optimally tracking targets[END_REF] considered the problem of tracking multiple mobile targets with multiple robots. They presented a general solution based on particle filtering in order to choose robot locations for the next time step that maximizes the quality of tracking. Frew [START_REF] Frew | Observer trajectory generation for target-motion estimation using monocular vision[END_REF] studied the problem of designing a robot trajectory, and not just the next robot location, in order to maximize the quality of tracking a single moving target. LaValle et al. [START_REF] Lavalle | Motion strategies for maintaining visibility of a moving target[END_REF] studied the problem of maintaining the visibility of a single target from a robot for the maximum time. Gans et al. [START_REF] Gans | Keeping multiple moving targets in the field of view of a mobile camera[END_REF] presented a controller that can keep up to three targets in one robot's field-of-view.

When the motion of the targets is fully known, the tracking problem can be formulated as a kinetic facility location problem. The goal of the stationary version is to place k facilities (robots) given the location of n sites (targets), so as to minimize the maximum distance between a facility and a site. For the kinetic version, Bespamyatnikh et al. [START_REF] Bespamyatnikh | Mobile facility location[END_REF] and Durocher [START_REF] Durocher | Geometric facility location under continuous motion[END_REF] presented approximation algorithms to control respectively one and two mobile facilities, when the trajectories for the sites are given. Recently, de Berg et al. [START_REF] Berg | Kinetic 2-centers in the black-box model[END_REF] presented improved approximation algorithms with two mobile facilities when only an upper bound on the velocities of the sites is available. However, the general problem of kinetic facility location with k facilities is open.

In the extreme case where no prior information of the targets is available, the multi-robot tracking problem can be formulated as a coverage problem [START_REF] Martinez | Motion coordination with distributed information[END_REF]. Schwager et al. [START_REF] Schwager | Eyes in the sky: Decentralized control for the deployment of robotic camera networks[END_REF] presented strategies to control the position and orientation of overhead cameras mounted on aerial robots in order to achieve equal visual coverage of the ground plane.

Unlike previous works, we study the trade-off between quality of tracking, and the number of targets tracked. We present an algorithm that chooses trajectories for each robot, instead of choosing just the next best location. This algorithm can be applied to the following two versions of the problem: tracking maximum number of targets, and maximizing the quality of tracking. We begin by formulating the problem and describing the sensing model.

III. PRELIMINARIES AND PROBLEM FORMULATION

Let k denote the number of robots, and n denote the total number of targets in the environment. The position of any robot or target is specified by their 3D coordinates x, y, z. The position of the i th robot at time τ is denoted by r i (τ ). Let z min be the minimum flying altitude. All robots have a camera that faces downwards. Let φ represent the field-ofview angle for the cameras. We assume that the robots can communicate amongst each other at all times.

Let t i (τ ) denote the position of the i th target. t i (τ ) is given by the position of a reference point that the robots can use to uniquely identify any target. For example, the reference point can be the centroid of a colored patch or a unique feature point on the object. All targets always move on the ground plane, i.e., z = 0 for all t i .

The reference point of any target t i in the field-of-view of a robot projects to some pixel in the image. A pixel can be backprojected to a ray in the world frame. In general, with no other information, it is not possible to solve for the target's location along this ray with a single camera measurement. However, since we assume that all targets move on the ground plane, we can solve for the coordinates of t i .

Ideally, we can exactly estimate t i given an image measurement, the camera pose, and the projection matrix. In practice, however, the following factors lead to an uncertain estimate of t i :

(1) The backprojection of camera pixels, which have quantized, integer coordinates, is not longer single ray but a pyramid (Figure 1(a)).

(2) Pixel measurements may be corrupted by noise. If the maximum noise is bounded by ∆p pixels, we backproject the set of pixels ±∆p around the measured pixel. The true target location is contained within the larger backprojection.

(3) The pose of the camera (or the robot) may not be accurately known. Typically, using exteroceptive sensors such as GPS and compass, we can bound the maximum uncertainty in estimating the robot pose. When the robot pose is known up to a bounded uncertain set, we can compute the backprojection for each pose within the set (Figure 1(b)).

In general, the quality of tracking under the three sources of errors, is a function of the relative distance and angle between the robot and the target, as seen in Figure 1(c). For a given true location of the target and an estimate of the robot pose, Figure 1(c) plots the maximum area of backprojection over all possible noisy measurements of the target, and all possible true robot poses.

While tracking, robots only have an estimate of the true target position. The uncertain estimate can be represented as a set of possible target locations on the ground plane. Given a motion model, the robots can propagate the set to obtain predicted target position, e.g., using particle filtering [START_REF] Thrun | Probabilistic robotics[END_REF]. The maximum area of backprojection can be computed for each predicted target position as shown in Figure 1(c).

The quality of tracking for a given target and robot pair can be defined as some measure of the areas of backprojection found for a predicted target position. Let q i (r j , τ ) denote the measure for target t i and robot r j at time τ . The quality of tracking t i at τ , is given by the best quality of tracking amongst all robots tracking t i , i.e., q i (τ ) = max j q i (r j , τ ). Finally, the total quality of tracking at τ is given by the sum of quality over all targets Q(τ ) = ∀i q i (τ ) over all targets. Alternatively, we may also consider the bottleneck quality over all targets Q(τ ) = min i q i (τ ).

IV. INFEASIBILITY OF TRACKING ALL TARGETS

In this section, we show the infeasibility of tracking all targets while maintaining any constant factor approximation of the optimal quality of tracking. We prove this by constructing an instance where the two goals, track all targets and maximize quality of tracking, conflict each other. We create a simple instance on a line where the quality of tracking is inversely proportional to the distance between the robot and the target: q i (r j , τ ) = 1/d(t i (τ ), r j (τ )) if t i is in the field-of-view of r j , and q i (r j , τ ) = 0 otherwise. The overall quality of tracking will be given by the bottleneck quality Q(τ ) = min i q i (τ ).

We use the instantaneous optimal quality of tracking, Q * (τ ), as the baseline for comparison. Q * (τ ) is the quality of tracking at τ , if one were to optimally place all the cameras at any location for any τ , regardless of their locations before τ . The placement of k cameras achieving Q * (τ ) may be significantly different from the placement achieving Q * (τǫ). There may or may not exist k continuous robot trajectories achieving Q * (τ ). Nevertheless, Q * (τ ) is an upper bound on the quality of tracking. This raises the question of whether we can at least maintain a constantfactor approximation of Q * (τ ) while tracking all targets. The theorem given next shows this is not possible, even when the motion of the targets is fully known.

Theorem 1 Let Q * (τ ) be the instantaneous optimal quality of tracking at time τ . Let the maximum speed of all targets be v. For any 0 < α ≤ 1 and β > 0, no algorithm can track all n > k targets with at least αQ * (τ ) quality for all τ with k ≥ 3 robots having a maximum speed of βv.

Proof: Consider Figure 2. We have k = 3 robots and n = 4 targets on a line. The distance between t 3 and t 4 is 0 at time 0. Targets t 1 , t 2 and t 3 remain stationary at all times, and t 4 moves with v = 1 to the right on the line. z min = 1 and φ = π/4 denote the minimum flying altitude and field-of-view angles (Section III). Fig. 2. At τ = 0, t 3 and t 4 are covered by the same robot to achieve Q * (0), where as for τ > d 12 , t 3 and t 4 are covered by separate robots.

If we have 4 targets and 3 robots, then there must exist a robot covering at least two targets at any given time. At τ = 0, we can verify that the optimal algorithm uses separate robots to cover t 1 and t 2 , and one robot to cover t 3 and t 4 (Figure 2). That is, Q * (0) = 1. Similarly, for any time τ > d 12 , optimal uses separate robots to cover t 3 and t 4 , and same the robot to cover t 1 and t 2 making

Q * (τ ) = √ 2
d12 . Thus, in any optimal algorithm, of the two robots covering t 1 and t 2 , one will switch to cover either t 3 or t 4 , after τ = d 12 . An approximation algorithm, on the other hand, does not necessarily have to make the same switch. Nevertheless, by setting d 12 appropriately, we will show that any approximation algorithm will be required to make the same switch at some time. By making d 23 sufficiently large, we will show that such a switch is infeasible with bounded velocity robots. The rest of the proof shows the existence of appropriate d 12 and d 23 values. This construction is similar to the one used by Durocher [START_REF] Durocher | Geometric facility location under continuous motion[END_REF] to prove the inapproximability of the kinetic k-center problem. For the case of aerial robots, however we show how to additionally take into account nonzero z min and φ values.

Let ALG be any algorithm that maintains a quality

Q(τ ) ≥ αQ * (τ ). If we set d 12 > √ 2
α , then ALG cannot use the same robot to cover t 1 and t 2 at time τ = 0. Else, Q(0) < α = αQ * (0) which violates the approximation guarantee. Hence, ALG uses separate robots to cover t 1 and t 2 at time 0.

Similarly, we can show that for any time τ > d12 α , ALG must use separate robots to cover t 3 and t 4 . Else Q(τ ) < √ 2 τ < αQ * (τ ) violating the approximation guarantee.

One of the two separate robots, say r, covering t 1 and t 2 initially, must cover either t 3 and t 4 at time τ > d12 α . In time τ , r must travel at least

d 23 -1 α -d12 √ 2α distance.
Here,

1 α and d12 √ 2α come from the condition that Q(0) ≥ α and Q(τ ) ≥ α √ 2 d12 . Consider a time τ = 2d12 α . At this time, r covers a maxi- mum distance of βτ = β 2d12 α . Set d 23 > β 2d12 α + 1 α + d12 √ 2α
. r cannot simultaneously cover at least one of t 1 or t 2 at time 0, and at least one of t 3 or t 4 at time τ , which is a contradiction. Hence, ALG cannot maintain an α approximation of Q * for all times.

The instance created in the proof above uses minimum flying altitude z min = 1 and camera field-of-view angle φ = π/4. We can create corresponding instances for any other values of these parameters. In light of Theorem 1, we drop the requirement that all targets must always be tracked. Instead we focus on the case when the robots are allowed to track a fraction of all targets.

V. 1/2 APPROXIMATION ALGORITHM In this section, we present the main algorithm to maximize the number of targets tracked, or maximize the quality of tracking. We divide the time into rounds of fixed duration. We consider the scenario where using measurements from previous rounds, the robots are able to predict the motion of the targets for the current round. For each robot, we create a set of m candidate trajectories that can be followed for the current round. For example, these trajectories can be generated using existing grid-based or sampling-based methods [START_REF] Green | Toward optimal sampling in the space of paths[END_REF]. Our goal is to choose a trajectory for each of the robots for the current round.

Figure 3 shows a simple instance with two robots, and three candidate trajectories each robot can follow. The camera footprint along two such trajectories as well as the set of targets covered by these trajectories are shown. Note that the trajectories need neither be restricted to any discretized grid, nor have uniform length or uniform speed.

Let R j (x) denote the set of targets predicted to be covered by x th trajectory followed by j th robot. We create a set system (X, R) where X is the set of all targets and R is a collection of all R j (x) sets. We group sets in R into k collections, one per robot. Each group contains m sets each. That is,

R = { R 1 (1), . . . , R 1 (m) candidate trajectories for r1 , . . . , R k (1), . . . , R k (m) candidate trajectories for r k } (1)
A valid assignment of trajectories can be represented by a map, σ : [1, . . . , k] → [1, . . . , m], indicating trajectory σ(j) (i.e., the set R j (σ(j))) is chosen for the j th robot. We can remove a target from the set R j (x) if it does not satisfy a given minimum quality of tracking requirement.

A. Maximizing Number of Targets

First consider the case of maximizing the number of targets tracked by k robots. This problem is a generalization of the maximum coverage problem [START_REF] Hochbaum | Analysis of the greedy approach in problems of maximum k-coverage[END_REF] stated as: choose k subsets to maximize the cardinality of the union of all subsets. In our case, we cannot arbitrarily pick k subsets since they must belong to distinct groups (i.e., the same robot cannot be assigned to two trajectories).

The maximum coverage problem, under group constraints, can be stated as: choose k subsets of R given by a map, Fig. 3. At the start of each round, we have a set of m candidate trajectories per robot. The trajectories may be non-uniform and of varying speeds. Using the predicted motion of the targets, we can determine which targets will be covered for a given trajectory and the corresponding quality of tracking.

σ : [1, . . . , k] → [1, . . . , m] such that the union of all subsets is maximized. The constraint that the same robot cannot be assigned to two trajectories is enforced by requiring the output be a map σ. This problem is known as the Maximum Group Coverage (MGC) problem. Chekuri and Kumar [START_REF] Chekuri | Maximum coverage problem with group budget constraints and applications[END_REF] proved that the greedy algorithm yields a 1/2 approximation for MGC. Their algorithm can directly be applied to track half the number of targets as an optimal algorithm. Our contribution is to extend the analysis to the weighted case.

B. Maximizing Quality of Tracking

For the case of maximizing the overall quality of tracking, we formulate a weighted version of MGC. Let q i (R j (x)) be the quality of tracking target t i with robot r j following the x th trajectory. q i (R j (x)) can represent the expected quality of tracking as described in Section III. The weight of any set R j (x) ∈ R is given by the sum of qualities of all targets tracked by R j (x). The objective is to maximize the sum of quality of tracking for all targets 1 .

The greedy algorithm for the unweighted MGC can be modified for the weighted setting (Algorithm 1). In each iteration, we choose a set R j (x) greedily that maximizes the total weight. We add R j (x) to the solution, and discard all other sets belonging to the same group, i.e., all other candidate trajectories for the same robot r j . This proceeds until we have chosen a trajectory for all robots.

Algorithm 1: Greedy Weighted MGC Algorithm

1 C ← ∅, I ← ∅ 2 for p = 1 to k do 3 Find R i (x) such that Q(R i (x) ∪ C) is greatest, and i ∈ I 4 σ(i) ← x 5 C ← C ∪ R i (x) 6 I ← I ∪ {i} 7 end 8 Return σ
Theorem 2 Algorithm 1 gives a (1/2ǫ) approximation for the weighted MGC problem for any ǫ > 0 in polynomial time. 1 The bottleneck version of maximizing the minimum quality of tracking over all targets cannot be applied since not all targets are tracked.

The analysis by Chekuri and Kumar [START_REF] Chekuri | Maximum coverage problem with group budget constraints and applications[END_REF] for the unweighted case can be modified for this weighted case. We present our full proof in the accompanying technical report [START_REF] Tokekar | Multi-target visual tracking with aerial robots[END_REF], for completeness.

We now evaluate the greedy algorithm through simulations and preliminary experiments.

VI. SIMULATIONS

In this section, we describe our implementation of the algorithm, and evaluate its performance through simulations. We carried out the simulations using the SwarmSimX simulation environment [START_REF] Lächele | SwarmSimX: Real-time simulation environment for multi-robot systems[END_REF]. SwarmSimX is a real-time multirobot simulator designed for modeling rigid-body dynamics in 3D environments. Models of the MikroKopter Quadrotor2 were used to simulate the motion of the robots.

For simulating the targets, we generated random trajectories as follows. Each target randomly chooses a speed and direction and moves along this direction for a random interval of time, drawn from a normal distribution. This class of trajectories is motivated by wildlife monitoring applications, where foraging animals have been found to follow such mobility models [START_REF] Benhamou | How many animals really do the levy walk?[END_REF]. The mean and standard deviation of the normal distribution were set to 10 s and 1 s, respectively in the simulations.

The target trajectories were restricted to 20 × 20 m square on the ground plane. The initial locations of all targets were chosen uniformly at random near the robot locations. A moving average filter of window length 5 running at 10 Hz was used to estimate the position and velocity of the observed targets for the next planning round. A measurement for a target was obtained only if it was contained within the fieldof-view of some robot.

For each robot, we created the following set of candidate trajectories: (a) stay in place, and (b) radially symmetric along 8 horizontal directions with a speed of 0.5 m/s. Thus, each robot could choose from a set of 9 trajectories in a round. Each round was set to a duration of 2 s. A trial consisted of 50 rounds.

Figures 4(a) and 4(b) show the effect of the number of robots and the maximum speed of the targets. As expected, the number of tracks and quality of tracking increases as the number of robots increase. Increase in the maximum speeds of the targets has the effect of spreading them further apart, which further reduces the number of targets that can be tracked. For these trials, the height of the robots was fixed to 3.5 m (i.e., the size of the camera footprint was fixed). Figure 5 shows the total number of targets tracked in one representative trial as a function of the time. Once the robots have lost track of a particular target, they do no receive any position information about that target. Thus, they cannot predict the future locations for a lost target, unless it appears again in the field-of-view of some robot.

For the simulations, we did not incorporate the uncertainty due to sensing. In the next section, we validate the uncertainty model and present results from a preliminary experiment using 4 aerial robots.

VII. EXPERIMENTS

In order to validate our sensing model and the algorithm, we performed trials on an indoor setup (Figure 6). The setup consisted of four quadrotors controlled using the TeleKyb framework [START_REF] Grabe | The TeleKyb framework for a modular and extendible ROS-based quadrotor control[END_REF]. All robots communicated directly with a central computer via a wireless XBee link. Each robot was fitted with a downward facing camera. The cameras streamed the live images wirelessly directly to the central computer. An indoor motion capture system was used for position feedback, while orientation is stabilized onboard.

A. Validating the Sensing Model

We first conducted trials to validate the sensing model presented in Section III. A robot was programmed to fly along a given trajectory at heights of 1 m and 1.5 m. The motion of the robot was smoothed, so as to ensure that the roll and pitch angles remained close to zero. Colored balls were placed on the ground (Figure 6). The pink and the yellow colored balls were fixed to motion capture markers to record their ground truth locations. All cameras were calibrated to obtain they camera parameters.

Figure 7 shows an image obtained using the on-board camera, along with the estimated and true locations of the balls. The backprojection area was computed considering ±5

B. Tracking Experiment

We implemented the greedy algorithm on the four robots. The controller on-board the robot was set to operate the robots smoothly in near-hovering mode at an average speed 0.5 m/s. Each round lasted for 3 seconds. The pink and yellow balls were moved manually (Figure 6). For this trial, (a)

2 -3 -2 -1 0 1 2 3 R 1 R 2 R 3 R 4 X(m) Y(m) t=110 secs -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -3 -2 -1 0 1 2 3 R 1 R 2 R 3 R 4 X(m) Y(m)
t=113 secs (a) At t = 110 s, R 3 chose the trajectory moving to the left to keep tracking the yellow target. the locations of the targets were obtained from the motion capture system. The robots used a moving average filter to predict the locations of the targets, based on previous measurements. A radius of √ 2 m was found empirically to correspond to the camera footprint when the robots operated at a height of 2.5 m. The robots had one of the four grid neighbors in the z = 2.5 m plane as candidate trajectories.

2 -3 -2 -1 0 1 2 3 R 1 R 2 R 3 R 4 X(m) Y(m) t=119 secs -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -3 -2 -1 0 1 2 3 R 1 R 2 R 3 R 4 X(m) Y(m)
Figure 8 shows the locations of the robots and the targets before and after two key rounds: at times 110 s and 119 s. The two rounds show events when the robots predicted that the target would move out of the coverage area in the next round. Hence, as an outcome of the greedy algorithm, the robots chose corresponding trajectories in order to continue to track the targets.

The sensing validation and tracking trials presented here demonstrate a proof-of-concept implementation of the components of our system. Our ongoing efforts are directed towards performing large scale experiments with this system.

VIII. CONCLUSION

In this paper, we studied a visual tracking problem in which a team of robots equipped with cameras are charged with tracking the locations of targets moving on the ground. We discussed the sources of uncertainty that affect the quality of estimating the locations of ground targets using overhead images. We showed the infeasibility of tracking all targets while maintaining the optimal quality of tracking, or any factor of the optimal quality, at all times. We then formulated the target tracking problem where the goal is to assign trajectories for each robot in order to maximize the quality of tracking. When we are given a set of candidate robot trajectories, we showed how the problem can be posed as a combinatorial optimization problem. A simple and easy-toimplement greedy algorithm applied to this problem yields a 1/2 approximation. Finally, we presented results from simulations and preliminary experiments validating the sensing model and demonstrating the feasibility of implementing the algorithm. Future work includes investigating the problem under inter-robot communication constraints, and conducting larger scale experimental validation.

Fig. 1 .

 1 Fig. 1. (a) Backprojection from a pixel yields a pyramid. (b) Uncertainty in target's estimate due to uncertain yaw angle of the robot. (c) Map showing the area of projection for the true target at [x, y, 0] (best viewed in color). The camera pose is estimated to have position [0, 0, 5] m and roll, pitch and yaw angles as 0 radians. Maximum image noise is ±5 pixels.

Fig. 4 .

 4 Fig. 4. (a) Number of targets covered out of 50 targets in the environment. (b) The average quality of tracking. The weight q i (R j (x) is computed as the inverse of the minimum distance between the target and the robot along R j (x).

Fig. 5 .

 5 Fig.5. The number of targets tracked in one trial. As the targets spread the total number of targets that can be tracked decreases. Once a target moves out of the field-of-view, the robots cannot predict their future locations.

Fig. 6 .

 6 Fig. 6. Experimental setup. Each robot is fitted with a downward facing wireless camera. All robots directly communicate with a central computer.

Fig. 7 .

 7 Fig. 7. Validating the sensing model. (a) On-board camera image. (b) The true target location (colored circles) in the global frame, and the estimated locations using the method described in Section III.

  t=123 secs (b) At t = 119 s, R 1 chose the trajectory moving to the right to keep tracking the pink target.

Fig. 8 .

 8 Fig. 8. Start (left figures) and end (right figures) of two rounds. Dashed trail shows the locations of the robots and targets in the preceding 5 secs.

Preprint version, final version at http://ieeexplore.ieee.org/ 2

Preprint version, final version at http://ieeexplore.ieee.org/ 3

http://mikrokopter.de Preprint version, final version at http://ieeexplore.ieee.org/

Preprint version, final version at http://ieeexplore.ieee.org/ 62014 IEEE IROS

IX. ACKNOWLEDGEMENTS

We thank the Max Planck Institute for Biological Cybernetics in Germany, for the experimental setup, and Simon Hartmann and Massimo Basile for their help during the experiments.

This material is based upon work supported by the National Science Foundation under Grant Nos. 1317788 and 1111638.