David Tschumperlé
email: david.tschumperle@greyc.ensicaen.fr

Luc Brun
email: luc.brun@greyc.ensicaen.fr

IMAGE DENOISING AND REGISTRATION BY PDE'S ON THE SPACE OF PATCHES

We define a simple transform able to map any multi-valued image into a space of patches, such that each existing image patch is mapped into a single high-dimensional point. We show that solving variational problems on this particular space is an elegant way of finding the natural patchbased counterparts of classical image processing techniques, such as the Tikhonov regularization and Lucas-Kanade registration methods. We end up with interesting variants of already known (non-variational) patch-based algorithms, namely the Non Local Means and Block Matching techniques. The interest of considering variational approaches on patch spaces is discussed and illustrated by comparison results with corresponding non-variational and non-patch methods.

INTRODUCTION

In the fields of Image Analysis, Processing and Synthesis, patch-based techniques generally meet with success. Defined as local square neighborhoods of image pixels, patches are very simple objects to work with, but they have the intrinsic ability to catch large-scale structures and textures present in natural images. Patches provide one of the simplest way to analyze, compare and copy textures, as soon as the considered patch sizes are higher than the so-called texel sizes (texture element, seen as the smallest significant unit of a texture). Moreover, patch-based methods are somehow intuitive : They indirectly reproduce the way humans are accomplishing some vision tasks by comparing semi-local image neighborhoods together. Proposed algorithms are often simple to implement, but they give surprisingly good results.

For instance, it has been long since patches have been used for solving the problem of estimating a displacement field between two images. This alignment problem can be classically (and partially) solved by the so-called Block Matching algorithm [START_REF] Haworth | Performance of Reference Block Updating Techniques when Tracking with the Block Matching Algorithm[END_REF][START_REF] Wiegand | Overview of the H264 Video Coding Standard[END_REF], which consists in comparing each neighborhood (patch) of one image with all of the other one, finding the best match for each pixel location. Unfortunately, the regularity of the resulting motion field is not ensured by such a crude method and regularized techniques [START_REF] Baker | Lucas-Kanade 20 Years On: A Unifying Framework[END_REF][START_REF] Lucas | An Iterative Image Registration Technique with an Application to Stereo Vision[END_REF] are sometimes needed. More recently, patch-based algorithms have been proposed to tackle the problem of synthesizing a texture similar to an input model [START_REF] Ashikhmin | Synthesizing Natural Textures[END_REF][START_REF] Efros | Texture Synthesis by Non-Parametric Sampling[END_REF][START_REF] Wei | Fast Texture Synthesis using Treestructured Vector Quantization[END_REF]. It has been shown that this quite complicated task can be very well achieved by a simple iterated copy-paste procedure of different patches coming from the model, found to fit the best with the local neighborhoods of the image to synthesize. Some variants of these techniques have been also applied for transferring textures from an image to another one [START_REF] Ashikhmin | Fast Texture Transfer[END_REF][START_REF] Hertzmann | Image Analogies[END_REF], and for textured inpainting [START_REF] Criminisi | Region Filling and Object Removal by Exemplar-based Image Inpainting[END_REF][START_REF] Hays | Scene Completion using Millions of Photographs[END_REF], i.e.the reconstruction of missing textured regions in an image. Despite their relative simplicity, one must admit that these algorithms give outstanding results. It is also worth to cite the set of recent patch-based denoising methods, initiated with the Non Local Means scheme [START_REF] Buades | A Non-Local Algorithm for Image Denoising[END_REF] and continued with various derivatives [START_REF] Azzabou | Variable bandwidth image denoising using image-based noise models[END_REF][START_REF] Boulanger | Space-time Adaptation for Patch-based Image Sequence Restoration[END_REF][START_REF] Brox | Iterated Nonlocal Means for Texture Restoration[END_REF][START_REF] Kervrann | Optimal Spatial Adaptation for Patch-based Image Denoising[END_REF]. Such methods are mainly based on an iterative weighted average of image patches. Here again, they have rapidly entered the hall of fame of image denoising techniques.

On the other side, variational methods for image processing [START_REF] Aubert | Mathematical Problems in Image Processing: PDE's and the Calculus of Variations[END_REF][START_REF] Bertalmío | Image inpainting[END_REF][START_REF] Sapiro | Geometric Partial Differential Equations and Image Analysis[END_REF] are known to be mathematically wellposed, flexible and competitive : They have the great ability to put complex a priori constraints (such as regularity constraints) on the obtained solutions. While regularity is desired for many applications (e.g.for motion estimation), it may sometimes avoid the reconstruction of textured solutions (e.g.for denoising/inpainting), textures being oscillating patterns by nature. Few attempts have been made to incorporate textured features in specific variational denoising/inpainting formulations by putting explicit patches [START_REF] Azzabou | Variable bandwidth image denoising using image-based noise models[END_REF][START_REF] Bertalmío | Simultaneous Structure and Texture Image Inpainting[END_REF][START_REF] Kindermann | Deblurring and Denoising of Images by Nonlocal Functionnals[END_REF][START_REF] Osher | Image Decomposition and Restoration using Total Variation Minimization and the H-1 Norm[END_REF] or various image transforms [START_REF] Fadili | Morphological Diversity and Sparse Image Denoising[END_REF] into the functionals to minimize, with few success. Mixing the performance of patches and the flexibility of variational methods is still a very exciting goal and an open challenge.

In this paper, we are looking toward a more general strategy to find patch-based counterparts of image processing techniques expressed as variational problems. Mainly, the idea lies on the construction of an alternate patch space on which the image to process/analyze is mapped. Thus, in this high-dimensional space, each existing image patch is represented by a single point. Variational formulations can be then defined and solved directly on the patch space, projecting back the obtained solution to the original image domain, if needed. So, instead of explicitly incorporating image patches into ad-hoc energy functionals, we rather consider the extension of existing classical (non-patch) functionals to a higher-dimensional space of patches. This is quite straightforward : Energy functionals are generally expressed with terms that can be easily extended for an arbitrary number of dimensions (e.g.gradients). As a result, we obtain algorithms which are the natural patch-based counterparts of "pointwise" variational techniques.

This paper is organized as follows : First, we define the reversible mapping of a multi-valued image into its patch space (section 2). Then, in this space, we consider the minimization of the Tikhonov regularization functional [START_REF] Bouhamidi | Sylvester Tikhonov-Regularization Methods in Image Restoration[END_REF][START_REF] Tikhonov | Regularization of Incorrectly Posed Problems[END_REF] (section 3). We show that the natural patch-based version of the resulting minimizing flow can be interpreted as a variant of the Non Local Means filter [START_REF] Buades | A Non-Local Algorithm for Image Denoising[END_REF]. In a second attempt, we tackle the problem of image alignment similarly, by minimizing an energy functional inspired by the Lucas-Kanade method [START_REF] Baker | Lucas-Kanade 20 Years On: A Unifying Framework[END_REF][START_REF] Lucas | An Iterative Image Registration Technique with an Application to Stereo Vision[END_REF]. We end up with an interesting variational version of the Block Matching algorithm, implemented by the evolution of semi-local nonlinear PDE's. Application results and discussions on possible future applications of this general variational framework on patch spaces conclude this paper (sections 5-6).

DEFINITION OF THE PATCH SPACE

Let us consider a 2D multi-valued image I : Ω → R n defined on a continuous domain Ω ⊂ R 2 . In this paper, we will mainly illustrate applications for n = 3, i.e.color images defined in the (R, G, B) color space. The i th component of I is a scalar image, denoted by I i : Ω → R, so that ∀(x, y) ∈ Ω, I (x,y) = I 1(x,y) , I 2(x,y) , . . . , I n(x,y) T .

More generally, the i th component of a vector X will be written X i and the restriction of X to a set of consecutive components i . . . j, as X i...j .

Patch definition :

The patch P I (x,y) located at (x, y) ∈ Ω on the image I is defined as the set of all image values belonging to a spatially discretized local p × p (square) neighborhood of I centered at (x, y). The spatial discretization step, related to the analysis scale of the considered patches, is assumed to be 1 in order to simplify notations, even though any step is possible. The size p is considered as odd, i.e.p = 2q + 1 (q ∈ N * ). Actually, a patch P I (x,y) can be ordered in a np 2 -dimensional vector as : P I (x,y) = I 1(x-q,y-q) , . . . , I 1(x+q,y+q) , I 2(x-q,y-q) , . . . , I n(x+q,y+q) T One may see P I (x,y) as the concatenation of the patch vectors P Ii (x,y) for all image channels I i , with P Ii (x,y) = I i(x-q,y-q) , . . . , I i(x+q,y+q)

T

An interesting mathematical study of the manifold formed by the set of all these patches have been initiated in [START_REF] Peyré | Image Processing with Non-local Spectral Bases[END_REF][START_REF] Peyré | Manifold Models for Signals and Images[END_REF]. This study is fascinating and share some ideas with the current paper.

Mapping to the patch space : We define the (np 2 + 2)dimensional patch space Γ = Ω × R np 2 . Each point p ∈ Γ is a high-dimensional vector whose coordinates may contain informations of any (x, y) location in Ω, as well as all values of any p × p patch P ∈ R np 2 . Obviously, in this patch space Γ, we want to highlight all the points of the form p = (x, y, P I (x,y) ), i.e.the points which precisely correspond to existing patches in I. In this context, we call p a located patch. We define then a function Ĩ in Γ such that Ĩ(p) is non zero only for these particular located patches of I :

Ĩ : Γ → R np 2 +1 , s.a. ∀p ∈ Γ, Ĩ(p) =      (P I (x,y) T , 1) T if p = (x, y, P I (x,y) ) T 0 elsewhere (1
) Thus, we call the application F such that Ĩ = F (I) a patch transform. Note that the value space of Ĩ has an extra component set to 1 for points at existing located image patches of I. This dimension can be compared with the one introduced when dealing with projective spaces : It plays a role of weighting action when inverting the patch transform F , i.e.retrieving back the multi-valued image I from Ĩ. Intuitively, it defines how much a located patch in Γ is meaningful, and by default, all existing patches of the original image I have the same importance.

In this paper, we want to show that it is worth to solve variational problems on Γ, expressed as the minimization of energy functionals E( Ĩ) rather than energies E(I) on the original image domain Ω. For this purpose, we have to point out the fact that Ĩ is a highly discontinuous function. To avoid derivation problems of the energies E, we will work on a continuous approximation Ĩǫ of Ĩ, where each original located patch p is not mapped into a single point, but into a normalized Gaussian function G ǫ with a variance ǫ close to 0. From a mathematical point of view, this defines the approximation as Ĩǫ = Ĩ * G ǫ which is C-infinite. In the sequels, with a slight abuse of notations and when no confusions are possible, we will denote Ĩǫ simply by Ĩ.

Back-projection on the image domain :

Due to the high dimensionality of Γ, there are of course no unique ways to compute the inverse transform of the patch-based representation Ĩ = F (I). We define a back-projection method based on two steps : First, we retrieve for every location (x, y) ∈ Ω, the most significant located patch P Ĩ sig(x,y) expressed in Ĩ. It is found to be the one with the maximum weight, i.e.P Ĩ sig(x,y) = Ĩ1...np 2 (x, y, P Ĩ max(x,y) ), with

P Ĩ max(x,y) = argmax q∈R np 2 Ĩnp 2 +1 (x, y, q T ) (2)
Note that if one perturbs only slightly the patch transform Ĩ of an image I, there are very good chances to find the most significant perturbed patches at the same locations as the original ones, i.e.P Ĩ max(x,y) = P I (x,y) , even though the pixel values of the modified patches may be different, i.e.P Ĩ sig(x,y) = P I (x,y) .

In a second step, the back-projected image I is reconstructed by combining these most significant patches together. Several strategies are possible : Here, we will use the simplest one, which simply consists in copying the center pixel of each significant patch P Ĩ sig(x,y) at its corresponding known location (x, y), normalizing it by its weight : ∀i ∈ [1, n], ∀(x, y) ∈ Ω,

I i(x,y) = Ĩip 2 + p 2 +1 2 (x, y, P Ĩ max(x,y) ) Ĩnp 2 +1 (x, y, P Ĩ max(x,y) ) (3)
More generally, we could have copied an entire sub-patch of P Ĩ sig at (x, y), while blending overlapped neighborhood patches according to their relative weights. This consideration on how to copy image patches appears frequently in the literature related to patch-based methods, and no "best" solutions have appeared yet. It is interesting to see that in our case, this choice is clearly delimited as a part of an inverse transform only. Now that the direct and inverse patch transforms are clearly defined, we study the application of some variational methods on the patch space Γ rather than on the original image domain Ω, for image denoising (section 3) and registration (section 4).

IMAGE DENOISING BY PATCH-BASED TIKHONOV REGULARIZATION

Suppose we have a multi-valued image I noisy : Ω → R n corrupted by some kind of noise. So will be its patch transform Ĩnoisy . We are looking for a patch-based minimizing flow able to regularize Ĩnoisy rather than process directly I noisy . For this purpose, we minimize the following energy E 1 , classically denoted as the Tikhonov regularization functional, which have been simply extended to the high-dimensional space Γ :

E 1 ( Ĩ) = Γ ∇ Ĩ(p) 2 dp (4)
where

∇ Ĩ(p) = np 2 +1 i=1 ∇ Ĩi(p)
2 is the habitual extension of the gradient norm for multi-valued datasets [START_REF] Zenzo | A Note on the Gradient of a Multi-image[END_REF]. Note that this multi-valued gradient includes also the gradient ∇ Ĩnp 2 +1 of the patch weights.

Minimizing flow :

The PDE flow that minimizes (4) is found by the derivation of E 1 ( Ĩ) using the Euler-Lagrange equations and by the expression of the corresponding gradient descent algorithm. It leads to the well known heat flow equation, which is in our case performed on the highdimensional patch space Γ :

   Ĩ[t=0] = Ĩnoisy ∂ Ĩi ∂t = ∆ Ĩi (5)
where ∆ is the Laplacian operator on Γ.

Similarly to denoising techniques based on classical diffusion PDE's [START_REF] Aubert | Mathematical Problems in Image Processing: PDE's and the Calculus of Variations[END_REF][START_REF] Barash | A fundamental relationship between bilateral filtering, adaptive smoothing, and the nonlinear diffusion equation[END_REF][START_REF] Barash | A common framework for Nonlinear Diffusion, Adaptive Smoothing, Bilateral Filtering and Mean Shift[END_REF][START_REF] Sapiro | Geometric Partial Differential Equations and Image Analysis[END_REF][START_REF] Sochen | Representation of Colored Images by Manifolds Embedded in Higher Dimensional Non-Euclidean Space[END_REF][START_REF] Tschumperlé | Vector-valued Image Regularization with PDE's: A Common Framework for Different Applications[END_REF], we are not particularly interested by the steady-state solution of (5), since it would roughly consist of a constant solution. We are rather looking for a solution of this multi-dimensional heat flow at a particular finite time t 1 . It has been already proven [START_REF] Koenderink | The Structures of Images[END_REF] that this solution is the convolution of the initial estimate Ĩnoisy with a normalized Gaussian kernel G σ with a stan-

dard deviation σ = √ 2 t 1 .
Here, this convolution has to be done on the high-dimensional patch space Γ :

Ĩ = Ĩnoisy * G σ , with ∀p ∈ Γ, G σ(p) = 1 (2πσ 2 ) np 2 +2 2 e -p 2 2σ 2 (6)
This is a quite straightforward way of defining the patchbased counterpart of the Tikhonov regularization process.

Interpretation in the image domain :

Interesting things arise when one tries to interpret the result of this patchbased heat flow [START_REF] Baker | Lucas-Kanade 20 Years On: A Unifying Framework[END_REF] in the original image domain Ω. In fact, the patch transform (1) tells us that Ĩnoisy vanishes almost everywhere in Γ, excepted on the points p = (x, y, P I noisy (x,y) ) T corresponding to the locations of the original image patches in I noisy . So, the convolution ( 6) can be simplified as : ∀(x, y, P) ∈ Γ, Ĩ(x,y,P) = Ω Ĩnoisy (p,q,P I noisy (p,q) ) G σ(p-x,q-y,P I noisy (p,q) -P) dp dq

We also notice that the locations of the most significant patches in Ĩ, as defined in (2), will be the same as the ones in Ĩnoisy , since the convolution (6) by a Gaussian kernel leaves the maxima of the patch weights Ĩnp 2 +1 unchanged [START_REF] Koenderink | The Structures of Images[END_REF]. So, the inverse patch transform (3) of Ĩ has a simple closed-form expression in Ω, and can be written directly from I noisy , considering that we have ∀(x, y) ∈ Ω, P Ĩ max(x,y) = P Using these relations together with (3), and after few calculus, we end up with :

∀(x, y) ∈ Ω, I (x,y) = Ω I noisy (p,q) w (x,y,p,q) dp dq Ω w (x,y,p,q) dp dq [START_REF] Barash | A common framework for Nonlinear Diffusion, Adaptive Smoothing, Bilateral Filtering and Mean Shift[END_REF] with w (x,y,p,q) = w s (x,y,p,q) w p (x,y,p,q) and w s (x,y,p,q) = 1 2πσ 2 e -(x-p) 2 +(y-q) 2 2σ 2 w p (x,y,p,q) = 1 (2πσ

Thus, each pixel value I (x,y) of the regularized image is the result of a weighted averaging of all noisy pixels I noisy (p,q) , the weight depending both on the spatial distance between points (x, y) and (p, q) (first term in ( 8)), as well as the similarity between corresponding patches centered on (x, y) and (p, q) (second term in ( 8)). Of course, this interpretation as a filtering process in Ω greatly simplifies the implementation step since it naturally suggests an algorithm in a lower-dimensional space Ω. But we have to keep in mind that it is actually the outcome of a classical Tikhonov minimizing flow on the patch space Γ. Note that the Gaussian functions in the spatial and patch-based terms of the weighing function w (x,y,p,q) (8) have the same standard deviation σ. Having different σ can be easily simulated by pre-multiplying the image I noisy by a ratio factor λ before processing, being then equivalent than having two different standard deviations σ patch = σ spatial /λ.

Link with other filtering methods :

The patch-based Tikhonov regularization method ( 7) is actually quite similar to the Non Local Means technique, as defined in [START_REF] Buades | A Non-Local Algorithm for Image Denoising[END_REF]. Differences are twofold : First, it is naturally defined as the solution of a minimizing flow acting on multi-valued images, while the original formulation has been expressed as an explicit nonlinear filter acting on scalar images. Second, the averaging weights w i,j of the Non Local Means algorithm are only based on the similarity between patches. Here, the weighting function ( 8) also considers the spatial distances between patches to average. Also, in the extreme case when the patch size p is reduced to 1 (so, a patch is just one point), the regularization method [START_REF] Barash | A common framework for Nonlinear Diffusion, Adaptive Smoothing, Bilateral Filtering and Mean Shift[END_REF] becomes the natural multi-valued version of another well know nonlinear image filter, namely the Bilateral Filtering method [START_REF] Barash | A fundamental relationship between bilateral filtering, adaptive smoothing, and the nonlinear diffusion equation[END_REF][START_REF] Barash | A common framework for Nonlinear Diffusion, Adaptive Smoothing, Bilateral Filtering and Mean Shift[END_REF][START_REF] Paris | A Fast Approximation of the Bilateral Filter using a Signal Processing Approach[END_REF][START_REF] Tomasi | Bilateral Filtering for Gray and Color images[END_REF]. These filters are known to be anisotropic in the image domain Ω. The benefit of minimizing the Tikhonov functional on the patch space Γ rather than on the image domain Ω is obvious in terms of regularization quality : It avoids the typical isotropic smoothing behavior of the Tikhonov regularization that usually over-smoothes the important image structures, such as the edges, corners and textures. A comparative figure of all these regularization algorithms is shown and commented in section 5.

IMAGE REGISTRATION BY PATCH-BASED VARIATIONAL METHOD

Let us now consider the problem of estimating a displacement field u : Ω → R 2 between two multi-valued images I t1 : Ω → R n (the reference image) and I t2 (the target image). This image alignment problem can be typically solved by a variational method, which aims at finding the vector field u that minimizes the following energy :

E 2 (u) = Ω α ∇u (p) 2 + D (p,p+u) 2 dp (9)
The user-defined parameter α ∈ R + imposes a regularity constraint on the estimated vector field u, if it chosen to be non zero. D (p,q) is a measure of the dissimilarity between image pixels I t1 (p) and I t2 (q) . Lot of different expressions for D have been already proposed in the literature [START_REF] Aubert | Mathematical Problems in Image Processing: PDE's and the Calculus of Variations[END_REF][START_REF] Baker | Lucas-Kanade 20 Years On: A Unifying Framework[END_REF][START_REF] Lucas | An Iterative Image Registration Technique with an Application to Stereo Vision[END_REF]. One of the most common choice makes the assumption that I t1 and I t2 are acquired under the same global illumination conditions (for instance, they are successive frames of a video sequence) and then, that the constraint of the Brightness Consistency holds. In this case, D can be reasonably chosen to be :

D 1(p,q) = I t1 σ(p) -I t2 σ(q) ( 10
)
where I t k σ = I t k * G σ are filtered versions of the images I t k , convolved by a normalized Gaussian kernel G σ . It allows the consideration of regularized reference and target images instead of possibly noisy ones. The choice D = D 1 in (9) leads to a variational extension of the Lucas-Kanade method [START_REF] Baker | Lucas-Kanade 20 Years On: A Unifying Framework[END_REF][START_REF] Lucas | An Iterative Image Registration Technique with an Application to Stereo Vision[END_REF] for multi-valued images.

Dissimilarity measure on the patch space : The dissimilarity measure D 1 can be naturally extended to deal with patches, by expressing it with the patch transforms Ĩt1 and Ĩt2 instead of the original images I t1 and I t2 . This is quite straightforward : We just replace the pointwise intensities by the most significant patches (2) present in the regularized versions of the patch transforms

Ĩt k σ = Ĩt k * G σ : D 2(p,q) = Ĩt1 σ(p,P Ĩt 1 max(p) ) -Ĩt2 σ(q,P Ĩt 2 max(q)
) (11)
Here, the smoothed patch transforms Ĩt k σ represent edgepreserving filtered versions of the reference and target images I t k , as mentioned in section 3. It means also that the localization of the most significant patches is known to be P Ĩt k σmax = P I t k . In the particular case where σ = 0, the dissimilarity function D 2 is the same as the one used for the classical Block Matching method. Moreover, when α = 0 (no regularization constraints are considered on u), the Block Matching is a global minimizer of [START_REF] Bertalmío | Simultaneous Structure and Texture Image Inpainting[END_REF].

Minimizing flow :

In a more general setting, the Euler-Lagrange derivation of (9) gives the set of coupled PDE's which locally minimizes the energy functional

E 2 : ∀j ∈ [1, 2], ∀x ∈ Ω,                      u [t=0] = 0 ∂u j(x) ∂t = α ∆u j + np 2 +1 i=1 Ĩt1 σi(x,P I t 1 (x) ) -Ĩt2 σi(x+u,P I t 2 (x+u) ) [∇G i ] j(x+u) (12) where
G i(x) = Ĩt2 σi(x,P I t 2 (x) )
. Here again, this minimizing flow ( 12) can be implemented directly in the image domain Ω, without having to explicitly compute and store the patch transforms Ĩt k σ . Indeed, it requires only the retrieval of the most significant patches in Ĩt k σ , which can be proven to be P Ĩt k σsig = P I t k regul , where I t k regul is the nonlinear filtered version of I t k by the Tikhonov regularization flow on the patch space [START_REF] Barash | A common framework for Nonlinear Diffusion, Adaptive Smoothing, Bilateral Filtering and Mean Shift[END_REF]. This patch-based registration PDE ( 12) is a local minimizer of a Block Matching-like objective function. Differences are twofold : First, it implicitly considers edgepreserving filtered versions of the reference and target images, instead of isotropically smoothed ones. But most of all, it is able to put an important a priori smoothness constraint on the estimated field u. This method combines then the significance of the patch description for the matching of image structures, while keeping the aptitude of the variational methods to impose useful constraints. These interesting properties are discussed and illustrated with a comparative figure in section 5.

APPLICATION RESULTS

We applied the different variational techniques on the patch space presented in this paper, on different color images considered in their original (R, G, B) color space.

Color image denoising : The Tikhonov regularization flow [START_REF] Barash | A common framework for Nonlinear Diffusion, Adaptive Smoothing, Bilateral Filtering and Mean Shift[END_REF] on the patch space Γ can be used to enhance degraded color images (or other multi-valued datasets). Fig. 1 compares it with the most connected algorithms, namely the Non Local Means, the Bilateral Filtering and the classical Tikhonov regularization performed on the image domain Ω. Synthetic white Gaussian noise (σ noise = 20) has been added to the original color image Barbara. For the honesty of the comparison, we have not applied the scalar versions of these related filters as defined in the papers [START_REF] Buades | A Non-Local Algorithm for Image Denoising[END_REF][START_REF] Tomasi | Bilateral Filtering for Gray and Color images[END_REF], but their multi-valued extensions. It gives indeed better denoising results than applying them channel by channel. The PSNR between the noise-free and restored images, as well as the parameters used for the experiments are displayed. For each method, the parameters have been manually tuned to optimize the obtained PSNR. As our proposed flow [START_REF] Barash | A common framework for Nonlinear Diffusion, Adaptive Smoothing, Bilateral Filtering and Mean Shift[END_REF] is actually very close in its final expression to the Non Local Means and Bilateral filters, the denoising results are of course very comparable in terms of quality. It seems to perform a little bit better on certain regions (some edges are sharper in the zoomed part). We believe it is due to the fact that our weighting function [START_REF] Bertalmío | Image inpainting[END_REF] considers both the similarity as well as the spatial distances between patches. The visual improvement is very subtle anyway. Additional results of our patch-based regularization technique are illustrated on Fig. 2. It is interesting to notice that in spite of its intrinsic isotropic nature, the Tikhonov flow on the patch space (5) clearly competes against sophisticated anisotropic diffusion PDE's for the smoothing of multi-valued images, such as the method proposed in [START_REF] Tschumperlé | Vector-valued Image Regularization with PDE's: A Common Framework for Different Applications[END_REF].

Color image registration :

The patch-based registration PDE's [START_REF] Brox | Iterated Nonlocal Means for Texture Restoration[END_REF] allow to estimate a smooth displacement field between two color images while considering a patch-based dissimilarity measure. Such a displacement field can be used then to generate a morphing video sequence, where new synthetic frames are added between the reference and target images (i.e.a re-timing process, where the reference frame corresponds to t = 0 and the target to t = 1). Fig. 3 shows a comparison of the re-timed results obtained by our proposed registration method [START_REF] Brox | Iterated Nonlocal Means for Texture Restoration[END_REF] and by the most related ones, i.e.the Block Matching algorithm and the pointwise version [START_REF] Bouhamidi | Sylvester Tikhonov-Regularization Methods in Image Restoration[END_REF] of the variational registration technique [START_REF] Bertalmío | Simultaneous Structure and Texture Image Inpainting[END_REF]. For each method, the estimated displacement field u is displayed as a colored image (first column). The mid-interpolated frame (i.e.reconstructed at t = 0.5) and a zoomed part of it are shown on the two last columns. This is indeed the hardest frame to reconstruct since it classically exhibits the largest ghosting artefacts due to motion estimation errors. To be able to handle large displacements, we considered a classical multiscale evolution scheme for the PDE's that minimize the pointwise and patch-based versions of the energy functional [START_REF] Bertalmío | Simultaneous Structure and Texture Image Inpainting[END_REF]. Without any surprise, the original Block Matching algorithm fails in reconstructing a smooth motion field, resulting in a unpleasant mid-interpolated frame reconstruction. The pointwise version of the variational registration technique ( 9) is quite competitive when the smoothness constraint is high enough (α = 0.1), but fails otherwise. Surprisingly, our proposed patch-based variational formulation [START_REF] Bertalmío | Simultaneous Structure and Texture Image Inpainting[END_REF] with the dissimilarity measure (11) outputs a quite good reconstructed frame, even when α = 0. Considering a patch-based dissimilarity measure seems to have a natural tendency to intrinsically regularize the registration problem. Anyway, best visual results are obtained with our variational patch-based registration method while considering regularization constraints on the displacement field u (last row).

CONCLUSIONS & PERSPECTIVES

In this paper, a reversible patch-based transform of multivalued images has been presented. More important, we illustrated the fact that considering variational methods in the corresponding patch space is a very nice way of designing the natural patch-based counterparts of classical image processing techniques. The obtained experimental results seems to confirm this impression. We strongly believe that this effort is a first step in reconciling both patch and variational worlds and that our proposed framework can be used to translate many other "pointwise" variational formulations into their natural patch-based counterparts. 12), and comparisons with Block Matching [START_REF] Haworth | Performance of Reference Block Updating Techniques when Tracking with the Block Matching Algorithm[END_REF][START_REF] Wiegand | Overview of the H264 Video Coding Standard[END_REF] and pointwise variational registration (9) methods.

Figure 1 .

 1 Figure 1. Illustration of the Tikhonov regularization on the patch space (7), and comparisons with most related filtering techniques.

Figure 2 .

 2 Figure 2. Additional results obtained by the Tikhonov regularization flow on the patch space (7).

Figure 3 .

 3 Figure3. Illustration of the variational image registration method on the patch space[START_REF] Brox | Iterated Nonlocal Means for Texture Restoration[END_REF], and comparisons with Block Matching[START_REF] Haworth | Performance of Reference Block Updating Techniques when Tracking with the Block Matching Algorithm[END_REF][START_REF] Wiegand | Overview of the H264 Video Coding Standard[END_REF] and pointwise variational registration (9) methods.


D Tschumperlé

Y Bentolila

J Martinot

J Fadili

FAST TIME-SPACE TRACKING OF SMOOTHLY MOVING FINE STRUCTURES IN IMAGE SEQUENCES

Keywords: Object tracking, Structure tensors, Integral lines

We address the problem of temporal tracking fine pointlike and filamentary structures exhibiting smooth motions in image sequences. By taking these specific restrictions into account, we put forward an original tracking method based on the search of integral lines in time-space structure tensor fields. The method is simple and very efficient regarding computation time and tracking precision, which allows a subpixel accuracy in both spatial and temporal domains. We suggest a numerical implementation of the algorithm which is based on a modified constrained Runge-Kutta scheme. Its performance and potential applicability are illustrated through two real cases : the tracking of internal linear features in composite materials observed with X-rays and the tracking of granular-shaped objects moving in a gas flow.

INTRODUCTION

Automatic tracking of moving objects in image sequences is a very challenging task and has a wide number of applications in digital image processing. Then this is not surprising that many different mathematical frameworks have been proposed so far in the literature to tackle the tracking problem, including now popularized Kalman filtering [START_REF] Jackson | Tracking deformable moving objects under severe occlusions[END_REF][START_REF] Peterfreund | Robust tracking of position and velocity with Kalman snakes[END_REF][START_REF] Terzopoulos | Active Vision, Tracking with Kalman Snakes[END_REF] and particle filtering [START_REF] Blake | Active Contours: The Application of Techniques from Graphics,Vision,Control Theory and Statis-tics to Visual Tracking of Shapes in Motion[END_REF][START_REF] Gordon | Novel approach to nonlinear/non-gaussian Bayesian state estimation[END_REF], among other techniques. In the instant case, we focus on a specific tracking application aiming at analyzing dynamic phenomena, which can be quantified by the measurement of smooth motions of fine structures in digitally acquired multispectral image sequences. We are mainly interested in the motion of pointwise or linear structures, such as fibers. This setting is usually encountered in a wide range of physical studies, including application fields in chemistry [START_REF] Pradel | Interaction between granular and gas flows: the pinning effect[END_REF], applied mechanics [START_REF] Tschumperlé | Wire Structure Pattern Extraction and Tracking from X-Ray Images for Riser Mechanisms[END_REF] or medical imaging [START_REF] Basser | Aldroubi In Vivo Fiber Tractography Using DT-MRI Data[END_REF]. Our application is not subject to real-time constraints, and the tracking algorithm proposed in this paper acts as a postprocessing method, having full access to forward and backward frame informations. However as the analyzed image sequences are generally composed of a huge amount of frames, we particularly paid attention to the cost of the algorithm in terms of computation time.

The paper is organized as follows : The preliminary section 2 describes first a simple but efficient spatial tracking technique which allows to retrieve filamentary structures respecting some smoothness properties in a single multispectral frame of the sequence. Then, we show that this idea can be extended to the spatio-temporal domain in order to track pointwise objects subject to smooth motions in the whole sequence (section 3). A modified 4th-order Runge-Kutta integration method [START_REF] Press | Runge-Kutta Method[END_REF] is thereby proposed to deal with the algorithm implementation (section 4). By the way, the combination of both spatial and temporal tracking methods can be used to track filamentary structures throughout an image sequence. We apply the suggested techniques on two real physical applications (section 5). The first one concerns the tracking of internal linear features in composite material observed on Xray images in order to a better understanding of mechanical forces in action. The second concerns the tracking of mixed granular objects moving in a gas flow for refining purposes.

SPATIAL TRACKING

We consider a multispectral 2D image I : Ω ⊂ R 2 → R n (typically, n = 1 for scalar-valued and n = 3 for color images), containing a filamentary object such as one fiber of interest we want to analyze. We denote by I i the ith vector channel of the image I. Given an initial user-defined seed point (x 0 , y 0 ) T belonging to the fiber, the goal of the tracking is to retrieve the function C : R → R 2 giving the spatial coordinates C (s) = (x (s) , y (s) ) T of all points composing the filamentary curve. The fiber is assumed to be smooth enough, and the tracking problem can be thereby formulated as finding the C 1 curve C such that :

   C (0) = (x 0 , y 0 ) T ∂C ∂s = u (C(s)) (1)
where s stands for the curvilinear abscissa of C and u (C(s)) ∈ S 1 is the unit tangent vector to the fiber, defined at least on each fiber point C (s) . The fiber is said to be tracked forward when s → +∞ and backward when s → -∞. A function integration step is obviously needed to retrieve the function C.

In order to determine the tangent vectors u (C(s)) , two natural assumptions are considered : First, the desired tracked object must be contrasted enough with respect to the image background. Second, we assume that this filamentary structure is composed of pixels having roughly the same vectorvalued intensity (color). For the scalar-valued case (n = 1), u should be typically related to the isophote direction ∇I T ∇I . In a more general manner, we suggest to handle the multispectral case (n ∈ N) and to gain robustness to noise by choosing u as the minimal eigenvector (unit vector associated with the smallest eigenvalue) of the so-called smoothed structure tensor G (C(s)) located as C (s) . G is defined as the field of 2 × 2 symmetric and positive-definite matrices :

G = i ∇I i ∇I T i * G σ (2)
where

∇I i = ∂Ii ∂x ∂Ii ∂y T
is the usual gradient vector of the channel I i , and G σ is a normalized 2D Gaussian kernel. The structure tensor is a natural generalization of the gradient for multispectral images [START_REF] Zenzo | A note on the gradient of a multi-image[END_REF] and has shown to have several good properties when used to retrieve multispectral isophote directions [START_REF] Tschumperlé | Fast Anisotropic Smoothing of Multi-Valued Images using Curvature-Preserving PDE[END_REF][START_REF] Weickert | Anisotropic Diffusion in Image Processing[END_REF]. G is defined on the whole domain Ω so as u, which can be thus seen as a vector field u : Ω → S 1 whose value on a fiber point C (s) gives the tangent vectors to C. Note that the smoothing parameter σ has a direct impact on the regularity of the eigenvector field u, and thereby implicitely controls the regularity of the tracked fiber C. In particular, higher smoothing may be used to get rid of fiber occlusion problems by reconstructing correct tangent vectors on occluded area. For simple reasons, we assumed here that u is correctly aligned. As an eigenvector represents an orientation without notion of direction, the computed field u can be in fact possibly discontinuous, independantly from the value of σ. This is anyway not an issue for the function integration step since we only need the orientation information, as the tracking direction is generally explicitely known (see Section 4). Moreover, it is interesting to notice that our filamentary structure tracking algorithm have in common ideas with some existing white matter tractography methods proposed in the field of diffusion-tensor MRI medical imaging [START_REF] Basser | Aldroubi In Vivo Fiber Tractography Using DT-MRI Data[END_REF]. The results of our filamentary structure tracking approach is illustrated on Fig. 1.

EXTENSION TO TEMPORAL TRACKING

Now, let us consider S : Ω × [0, T ] → R n as the multispectral image sequence containing pointwise objects in motion we want to track. S is a volume W × H × T of stacked multispectral image frames. We limit ourselves to pointwise objects tracking since we are able now to detect filamentary structures starting from a single point (see section 2). Similarly to the spatial case, the tracking problem is formulated as follows : considering an user-defined input point P 0 = (x0, y0, t0) in one frame t0 of the sequence S, find the positions P (s) = (x (s) , y (s) , t (s) ) of this point in all frames, i.e at least when t (s) ∈ N. We assume that the tracked pointwise object is contrasted enough and exhibits a smooth motion throughout the sequence. Then, the trajectory P : [0, T ] → R 3 can be computed as the integration of a tangent vector field

w : Ω × [0, T ] → R 3 :    P (0) = (x 0 , y 0 , t 0 ) T ∂P ∂s = w (s)
This kind of tracking can be actually seen as a spatio-temporal tractography. The main challenge here consists in the choice of the tangent vector field w. Unfortunaly, the temporal extension of the tracking method is not as simple as choosing w to be the minimal eigenvector of the structure tensor field (2) of S (which is now a 3 × 3 matrix). Intuitively, it would correspond to compute the direction of smallest intensity variations in the image sequence viewed as a volumetric image. But nothing ensures that this 3D direction would have a non null temporal component. In fact, most of the smallest variations are often encountered in the same frame as the current tracked point, since objects we want to track are not perfectly pointwise but are usually composed of regions with pixels having similar intensities. This is for instance the case for fibers in Fig. 1 and balls in Fig. 2. As a result, the minimal eigenvector of the 3 × 3 structure tensor is generally not representative of the object motion only, but also of its spatial structure. Thus, mixing variations of all spatial and temporal axes (x, y, t) in the same structure tensor is definitely not a good idea.

The solution relies on the simultaneous but separated tracking of the considered point in the (x, t) and (y, t) spaces instead of considering the whole 3D domain (x, y, t). We propose to deduce the spatio-temporal tangent vector w by analyzing the 2 × 2 smoothed structure tensor fields G xt , G yt both defined on the domain Ω × [0, T ] and computed as

G xt = ∇S ixt ∇S T ixt * G σ and G yt = ∇S iyt ∇S T iyt * G σ where ∇S ixt = ∂Si ∂x ∂Si ∂t T , ∇S iyt = ∂Si ∂y ∂Si ∂t T
are the two spatio-temporal gradients respectively in (x, t) and (y, t) spaces. For each point (x, y, t) of the sequence, the minimal eigenvectors u xt = (u xt , v xt ) T and u yt = (u yt , v yt ) T of G xt and G yt give coherent measures on the spatio-temporal variation directions of S. Mixing the directions u xt and u yt allows to retrieve a good approximation of the point motion, and ensures that the estimated motion tangent vector w goes through time :

w =      (u xt , uyt|vxt| |vyt| , v xt ) T if |v xt | < |v yt | and v xt = 0 ( uxt|vyt| |vxt| , u yt , v yt ) T if |v xt | > |v yt | and v yt = 0
Intuitively, the tracking procedure reads as follows. For each curve point (x (s) , y (s) , t (s) ), we get two motion vectors u xt and u yt in the (x, t) and (y, t) spaces, characteristic of how important is the estimated temporal motion along x and y. The final motion vector w is then chosen mainly from the one that varies the less, and is extended to the 3D domain using the data provided by the other. It appears that this relatively simple motion vector computation works surprisingly well in order to estimate the coherent motion appearing in the image sequence.

IMPLEMENTATION CONSIDERATIONS

The main issues concerning the numerical implementation of the tracking algorithm are twofold.

Computation of the 2 × 2 structure tensors : A classical way of estimating the structure tensor field consists in computing a smoothed version of

G = n i=1 ∇I c i ∇I c T i where ∇I c i = 1 2 (I i(x+1,y) -I i(x-1,y) , I i(x,y+1) -I i(x,y-1) )
T is the classical central finite difference scheme of the image gradient. Here, we would rather suggest a slightly modified scheme which advantagely uses the more precise forward and backward finite difference schemes, while remaining symmetric :

G = n i=1 ∇I f i ∇I f T i + ∇I f i ∇I b T i + ∇I b i ∇I f T i + ∇I b i ∇I b T i 4
where ∇I f i and ∇I b i correspond to the classical forward and backward differences :

∇I f i = I i(x+1,y) -I i(x,y) I i(x,y+1) -I i(x,y)
and

∇I b i = I i(x,y) -I i(x-1,y) I i(x,y) -I i(x,y-1)
The same kind of formulae is used to compute the spatiotemporal structure tensors G xt and G yt .

Curve integration :

The problem consists in numerically integrating the equation (1). A classical Euler scheme would consist in iterating :

C (s+ds) = C (s) + ds w (C(s))
with ds, an user-defined integration step close to zero. A more precise integration algorithm has been obtained here, using the so-called 4th-order Runge Kutta scheme [START_REF] Press | Runge-Kutta Method[END_REF]. It consists in introducing several interpolations of the motion vector w for one iteration :

k 1 = ds u(C (s) ), k 2 = ds u(C (s) + k1 2 ), k 3 = ds u(C (s) + k2 2 ), k 4 = ds u(C (s) + k3 2 )
Then, the next point C (s+ds) is obtained by

C (s+ds) = C (s) + k1 6 + k2 3 + k3 3 + k4 6
Linear interpolation has been used successfully to estimate the value of the vector w on inter-pixel coordinates.

As noticed in section 2, only the orientation of w is relevant, not its direction. To get rid of this direction, we have introduced a realignement step before each tracking iteration, such that the estimated vectors w in the 3 × 3 neighborhood of the current point C (ds) are aligned with the one used for the previous iteration. It ensures that the tracking is done in a coherent direction without half-turns.

APPLICATIONS AND RESULTS

Some results of our proposed spatial and temporal tracking algorithm is illustrated with two different applications, one for mechanical purposes with a sequence containing fibers, and the other for refining purposes, with a sequence of moving pointwise granular-shaped objects.

Fiber Tracking : In [START_REF] Tschumperlé | Wire Structure Pattern Extraction and Tracking from X-Ray Images for Riser Mechanisms[END_REF], we have already described a complete image processing pipeline in order to separate and display internal structures of a composite material submitted to external stresses. Fibers composing this material are moving and must be tracked to understand their mechanical behavior. Image sequences of these moving fibers are acquired thanks to a specific X-ray video camera, then filtered using techniques presented in [START_REF] Tschumperlé | Wire Structure Pattern Extraction and Tracking from X-Ray Images for Riser Mechanisms[END_REF], including a denoising step [START_REF] Tschumperlé | Fast Anisotropic Smoothing of Multi-Valued Images using Curvature-Preserving PDE[END_REF], a waveletbased image contrast enhancement step [START_REF] Starck | Gray and Color Image Contrast Enhancement by the Curvelet Transform[END_REF], followed by a Morphological Component Analysis (MCA) [START_REF] Starck | Image decomposition via the combination of sparse representations and variational approach[END_REF], enabling to distinguish between fiber components of the image. The spatial and temporal tracking algorithms described in section 2 and 3 have been succesfully applied to extract the shape of the fiber (Fig. 1).

Granular-shape tracking : This example illustrates a refinery process, which consists in injecting a gas across a catalyst moving bed (complete description of the experiment is detailed in [START_REF] Pradel | Interaction between granular and gas flows: the pinning effect[END_REF]). It involves a complex medium : a mix of granular material with a flowing gas inside. We have been interested in a phenomena which lies inside the refinery reactor when a horizontal flow goes through a moving granular bed. Thanks to image analysis applied to experimental equipment, Dynamic features are extracted to characterize the cavitation phenomena; it is based upon the kinematics of the grains. Here, our algorithm enables the tracking of a granular-shape object through a gas flow. The results are illustrated in the image sequence in Fig. 2.

CONCLUSION & PERSPECTIVES

We have proposed an original spatio-temporal tracking algorithm which is very simple to implement. It does not need expensive computation time and gives tracking results with a sub-pixel precision. Detailled numerical schemes have been proposed in order to get rid of the most common implementation issues. We have illustrated its interest with two various applications in both fields of chemistry and applied mechanics.

Fig. 1 .

 1 Fig. 1. Tracking of one fiber in a X-Ray sequence.

Fig. 2 .

 2 Fig. 2. Tracking of one granular-shaped ball in a color sequence.


D Tschumperlé
email: david.tschumperle@greyc.ensicaen.fr

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

Greyc Image

Y Bentolila
email: yohan.bentolila@ifp.fr

Wire Structure Pattern Extraction and Tracking From X-Ray Images of composite Mechanisms

This paper introduces a complete pipeline of image processing methods in order to analyze and track the internal structures of a composite material. As a first step, input Xray images are denoised, enhanced, separated into different morphological components, and geometrically filtered in order to isolate the interesting fiber inside the composite material. This requires the design of specific algorithms preserving very thin image details while being able to remove undesired image regions that may be sometimes large. For this purpose, we use state-of the-art techniques based on recent achievements in diffusion PDE's and modern harmonic analysis tools. Then, the shapes of the remaining fibered composite are individually analyzed by the use of a tensor-based tracking algorithm. We illustrate how this set of techniques allows the dynamic analysis of the composite structure when submitted to external mechanical loads.

Introduction and Context

Microstructure study of Composite material under stresses involves many complex mechanical phenomena and raises a strong interest in the off-shore industry [START_REF] Odru | Composite Riser And Export Line Systems For Deep Offshore Applications[END_REF][START_REF] Toussaint | Composite Riser for Ultra deep Offshore: the CORFU Project[END_REF]. Composite structure used as risers are typically made of many layers. Usually composite material is made of polymer reinforced with glass or carbon fibers. Due to the use of X-ray technology, the composite fiber material has been replaced by steel in order to visualize the fiber within the composite. Fiber layers are set with opposite lay angles and are basically used to sustain tensile loads. A way of understanding the mechanical phenomena applying on a composite part submitted to external stresses is to monitor the behavior of its various layers using an X-ray camera. For this purpose, devoted mechanical testing equipment has been developed between IFP and CEA 1and allows to stress the composite periodically, while recording the behavior of fiber layers using a numerical Xray video recorder (Fig. 1). The analysis of acquired X-ray images points out several properties which yield to develop appropriate image processing methods to extract patterns of interest. First, the used X-ray acquisition system adds non-negligible Poisson noise to the data. Second, the fiber patterns of interest have a quite low contrast while being at the same time only few pixels large. Third, black numbers corresponding to reference elements on the composite structure, hide fiber structures that exist behind them. Fourth, vertical structures are clearly visible on the image. They actually correspond to the underlying metallic frame of the experimental device and thus, it is caught by the X-ray acquisition. Last but not least, front and back layers appear as mixed in the acquired image. This is due to the fact that the X-ray system acquires the 2D projection of the composite sample (which is initially a 3D object). This yields to a lot of cross line structures in two mainly opposite directions. With such conditions, extracting the shapes of fiber structures with accuracy in this harsh scene is a challenging goal that requires advanced image processing tools. In this paper, we propose a set of image enhancing and analysis techniques that are sequentially applied in order to track such fiber patterns in good conditions. All these methods will be illustrated on a real dataset along the article.

Proposed image processing pipeline

Overview of the image processing chain

As pointed out in the introduction, the data of wire structures are not directly accessible, since a lot of degradations, artefacts and information mixing are present in the acquired X-ray images. We first tackle the denoising problem (section 2.2) by the use of a specific curvature-preserving regularization PDE, as proposed in [14]. As we want to enhance and isolate fiber structures in our images, noise removal has to be the first process to perform. It simultaneously avoids the enhancement of undesirable image structures and reduces the need of regularizing terms in image processing techniques used afterwards. In a second step, a wavelet-based image contrast enhancement (section 2.3) is applied in order to put forward fiber structures and prepare the image data for a source separation step using the so-called Morphological Component Analysis (MCA) (see section 2.4). This morphological decomposition of the image allows to clearly separate the image into two morphologically distinct components: one piece-wise smooth and the other with linear-curvilinear structures. An additional orientation filtering scheme is finally applied in order to remove vertical stripes due to the main metallic frame of the experimental device, and to separate the front and back fiber layers of the considered object (section 2.5). As a result of this image processing chain, two images are obtained corresponding to the two wire structure directions. In a last step, a fiber tracking process is performed to extract the wire pattern trajectory (section 2.6).

Image denoising

Let I : Ω → R be one of the acquired X-ray image of the considered composite. In order to remove the high amount of noise present in I, we used a technique based on regularization PDE's (partial differential equations), which are well known to be able to remove noise while preserving important discontinuities of the image data (by the mean of anisotropic filtering) [START_REF] Aubert | Mathematical Problems in Image Processing: PDE's and the Calculus of Variations[END_REF][START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF]9,17]. Here, important image data consist mainly in fiber structures on the composite sample. There are usually structures that are only few pixel large. In order to take this particularity into account, we use the framework of curvature-preserving PDE's, as proposed in [14]. We applied the following equation on the input Xray images :

∂I ∂t = trace(TH) + 2 π ∇I T π α=0 J √ Taα √ Ta α dα (1)
where a α = (cos α sin α) T , H and J √ Taα are respectively the Hessian of I and the Jacobian matrices of √ Ta α (matrices of the first and second spatial derivatives), while T is the field of 2 × 2 diffusion tensors :

T = 1 (λ++λ-) p + θ + θ T + + 1 (λ++λ-) p -θ -θ T -
Here, λ + , λ -and θ + , θ -respectively denotes the eigenvalues and eigenvectors of the smoothed structure tensor field

G = (∇I∇I T ) * G σ [17].
User-defined parameters of the algorithm p + and p -defines the local smoothing strength. This recently proposed PDE-based method has the advantage of preserving wire curved structures in the images, thanks to the constrained term J √ Taα √ Ta α added to the second order elliptic term trace(TH) that basically performs a diffusion along the tensor field T [START_REF] Cabral | Imaging vector fields using line integral convolution[END_REF]14]. Moreover, this method can be implemented with a smart and fast numerical scheme, based on Line Integral Convolutions (as defined in [3]). It allows to perform image denoising with a sub-pixel accuracy, which is particularly well suited for our aim of image denoising with preservation of thin fiber structures. A result of our image denoising method is illustrated on Fig. 2.

Contrast enhancement

The contrast enhancement technique that we propose is based on multiscale transforms, and more specifically on the wavelet transform (WT). Basically, the WT can be seen as a sophisticated multi-scale differential operator that will capture the edges of the image. Then, a non-linear enhancement function is applied coefficient-wise to each wavelet sub-band, and an inverse wavelet transform is finally applied to get the enhanced image. Several approaches have been developed in the literature for multiscale image enhancement (see e.g. [START_REF] Starck | Gray and Color Image Contrast Enhancement by the Curvelet Transform[END_REF]6] and references therein). Here we use a non-linear global adaptive gain (GAG) function that will enhance the significant high WT coefficients and reduce the non-significant small ones (mainly due to noise). By global, we mean that the GAG function will have the same expression independently of the wavelet coefficient location. Towards this goal, the following GAG function was designed:

f (d j,k ) = d j max a sigm(c(d j,k -b j )) -sigm(c(d j,k + b j )) (2) sigm(x) is the sigmoid function (1 + e x ) -1 . d j,k = d j,k /d j
max , where d j,k is the wavelet detail coefficient at scale j and location k, and d j max is the maximum wavelet coefficient in magnitude at scale j. c is the gain parameter (i.e. higher value of c yields more details enhancement). a and b j are defined as follows:

a = (sigm(c(1 -b j )) -sigm(c(1 + b j ))) -1 , b j = T d j max
b j can be thought of as a "threshold" below which the wavelet coefficients are supposed to be negligible. Therefore, a possible meaningful value would correspond to T = σ 2 log N 2 , which is best known as the universal threshold of Donoho and Johnstone. More details can be found in [6].

The results are depicted in Fig. 3, where the fine detail structures have been enhanced. Some faint structures that were hard to distinguish before this step are now clearly visible.

Source separation

Morphological Component Analysis (MCA) is a recent novel approach in modern harmonic analysis proposed by Starck et al. [START_REF] Starck | Image decomposition via the combination of sparse representations and variational approach[END_REF]. Its primary goal is image separation into several semantic parts (sources), where each component is sparsely represented in a given image transform. Suppose that an image I defined on a compact is a linear superposition of several parts {S i } i=1,...,ns , where each part has a parsimonious decomposition in a transform Φ i . The union of all these transforms is usually called a (possibly hyper-) redundant dictionary. For example, if the image is composed of a piece-wise smooth part and a textured part (the so-called cartoon+texture model [16]), then two very competitive candidate transforms in the dictionary would be the wavelet transform (for the piece-wise smooth part) and the local DCT (for the locally oscillating texture). Using both sparsity and variational mechanisms, the MCA approach attempts to solve the following optimization problem:

min {αi}i=1,...,n s J(α i ) = 1 2 X - ns i=1 Φ i α i 2 2 + λ ns i=1 α i 1 +µ S j T V
Here, S j = Φ j α j is known to be the piecewise regular component, and . T V is the Total Variation norm. Penalizing with TV, will force the component S j to be closer to a piecewise smooth image and, thus, can be helpful to support the separation process. In this formulation, beside the data fidelity term, the l 1 norm penalty on the transform coefficients promotes the sparsity of the representations.

Optimization issues

In [START_REF] Starck | Image decomposition via the combination of sparse representations and variational approach[END_REF], the authors proposed a modification of the Block Coordinate Relaxation (BCR) algorithm to solve the above optimization problem. The interested reader may refer to [START_REF] Starck | Image decomposition via the combination of sparse representations and variational approach[END_REF]10] for further details.

Candidate dictionaries

Our goal is to separate the X-ray denoised and enhanced into two components: one containing the curvilinear fiber structures and the other with the lead markers appearing as black numbers in the X-ray image. These two components will be denoted hereafter S c and S w . Therefore, natural candidate transforms in the dictionary would be the curvelet transform [4] for S c , which is very efficient to capture curvilinear structures, and the undecimated wavelet transform for the numbers and the piecewise constant background component S w . The TV penalty will be added to direct the image S w to fit the piecewise smooth model.

Orientation filtering

As stated in the introductory part, the goal of this step is twofold: (a) separate the wire structures according to their orientation in the component S c obtained from the MCA step, (b) remove all linear structures whose main orientation is either horizontal or vertical. This amounts to a "geometrical" filtering of the component S c . We then use the curvelet coefficients to select those coefficients that exhibit the desired geometrical (orientation) properties, and kill the others. The inverse curvelet transform is applied only to the retained coefficients to get the orientation-filtered image(s). More specifically, the component S c is efficiently represented by the curvelet transform, and can be written as:

S c (x, y) = γ α c γ φ γ (x, y)
where φ γ are the curvelet transform atoms indexed by their scale j, orientation θ and position (k, l) (that is γ = (j, θ, k, l)). An orientation-filtered version of the image S c , where all orientations are kept except θ * can be written:

S * c (x, y) = γ∈Γ * α c γ φ γ (x, y)
where Γ * = {(j, θ, k, l) s.t. θ = θ * }, i.e. Γ * is the set of indices γ containing all orientations but θ * . For example, in order to remove the horizontal and vertical structures, θ * takes the values 0, π and ±π/2. As far as the fiber structures are concerned, they have global orientations at π/4 and 3π/4. A natural choice for the angle θ * is obviously ±π/4 and ±3π/4. However, because of classical discretization issues induced on a Cartesian grid, and also because there can be some local fluctuations on the fiber orientation, some tolerance must be allowed on θ * to effectively remove all undesirable oriented structures. Furthermore, the number of orientations in the directional part of the curvelet transform is scale-dependent (see [4] for details). Thus, the angle tolerance on θ * must also depend on the scale j. Formally, this corresponds to the set Γ * = {(j, θ, k, l) s.t. θ = θ * ±∆θ * j }, where ∆θ * j is the scale-dependent angle tolerance. In our application, we set ∆θ j = 2π/n j θ where n j θ is the number of angular directions at scale j. Results obtained after orientation filtering are shown in Fig. 5. One important consequence of this step is that the filtered image is free from fiber structures crossings and undesired horizontal and vertical structures. This has an important impact on wire structures tracking.

Wire Structure Tracking

Having isolated the image data corresponding to interesting fiber structures on the composite, we are now able to retrieve more easily the fibers shape, since a fiber is now represented by an isovalued curve of the image pixels. Starting from an user-defined point belonging to a fiber of interest, we get its global shape by tracking the curve of constant image intensities (the color of the seed point) in both forward and backward directions. This is done by computing the integral curve of the vector field w : R → R 2 defined as the field of minimal eigenvectors of the smoothed structure tensor G = (∇I∇I T ) * G σ . This curve integration is performed through a 2nd order Runge-Kutta scheme [8].

Using this integration method allows us to find the path of minimal variations (i.e constant intensity) starting from a point , while introducing an interesting smoothness constraint on the tracked curve by the parameter σ which is the standard deviation of the structure tensor gaussian smoothing. This regularity prior is justified by the nature of the tracked object. Results of fiber tracking in one frame (back and front layers) are illustrated on Fig. 6.

Figure 6. Result of fiber tracking on front and back layers of the composite (green line).

Our scheme is actually similar to the one used in recent works on fiber tracking in 3D volumetric images of diffusion tensor MRI [START_REF] Tschumperle | Diffusion tensor regularization with constraints preservation[END_REF].

Figure 1 .

 1 Figure 1. One of the acquired X-ray image of a composite sample.

Figure 2 .

 2 Figure 2. PDE-based denoising of X-ray images (details) : Acquired images (left), PDE-regularized images (right).

Figure 3 .

 3 Figure 3. Effect of multi-scale contrast enhancement applied to a PDE-denoised image: top: enhanced image, bottom: zoom of an image region before (left) and after (right) enhancement.

Fig. 4

 4 Fig.4provides an illustration of the separation provided by the MCA for the above dictionaries after 30 iterations. The top row depicts the PDE-denoised and multiscale enhanced image, the middle and bottom rows give the two components S c and S w . The MCA performed impressively well. One can clearly see how it managed to get rid of the numbered markers while preserving the curvilinear fibers structure, and even reconstructing the unobserved fibers parts that were behind the markers.

Figure 4 .

 4 Figure 4. Results of the MCA separation step. Top row: PDEdenoised and multiscale enhanced image. Middle and bottom rows: components Sc and Sw provided by the MCA.

  (a) Wire structures after MCA (Sc image). (b) Image of (a) after horizontal and vertical lines were filtered out. (c) Image of wire structures in the NW and SE directions. (d) Image of wire structures in the NE and SW directions.

  (e) Zoom of (a),(b),(c),(d).

Figure 5 .

 5 Figure 5. Results of orientation filtering applied to the Sc MCA component.

French Atomic Energy Comission (CEA).

Acknowledgement

The development of this image processing method has been performed within an industrial project sponsored by Technip. The authors would like to express their gratitude to Technip research and development department for facilities access and data supports.

Conclusion & Perspectives

In this paper, we combined elements from PDE-based image processing and modern harmonic analysis to provide an effective solution to a challenging and complex image pro-cessing problem: fiber structures tracking on a composite sample observed at the output of an X-ray image acquisition system. Because of the low quality X-ray images, and many other complicating factors (e.g. projection of 3D structure onto a 2D plane), direct tracking on the original image is not possible. We then proposed to enhance the quality of the image and isolate the fiber structures prior to the fiber tracking step. Towards this goal, a curvature-preserving PDE-based approach was used for image denoising, while harmonic analysis with variational tools were proposed to enhance and separate the fiber structures. Finally, a tensorbased filamentary structure tracking algorithm was applied on the enhanced and simplified image. The obtained results are very precise and the complete pipeline of proposed image processing methods was successfully validated on a large database of composite X-ray images. Nonetheless, at this time, the tracking algorithm is still at a preliminary stage and deserves a much deeper work. Our ongoing work is directed towards this aspect.


Saleh Mousa Alzahrani
email: salzahrani@uqu.edu.sa

Xavier Antoine
email: xavier.antoine@univ-lorraine.fr

Chokri Chniti
email: cachniti@uqu.edu.sa

Computation of multiple scattering problems by circular cylinders

Keywords: Multiple scattering, integral equations, numerical simulation AMS classification. 76Q05, 78A45, 31A10, 65M70

The aim of this paper is to develop some well-adapted formulations for numerically solving the scattering problem of an incident time-harmonic wave by a cluster of circular cylinders. These formulations are illustrated by numerical simulations for different physical situations involving various kinds of boundary conditions.

Introduction

Multiple scattering by clusters of obstacles is a complex physical and computational challenge. It indeed results from the nontrivial scattering interactions between the single scatterers that compose the structure, yielding then to complex patterns of the wave field that are difficult to predict. Related to this specific behavior emerge some new physical properties of waves that are used to create some technological devices or to explain some observations in the nature. From the perspective of numerical methods, quite a lot of recent computational approaches were developed over the years. We refer e.g. to [START_REF] Acosta | On-surface radiation condition for multiple scattering of waves[END_REF][START_REF] Acosta | Coupling of Dirichlet-to-Neumann boundary condition and finite difference methods in curvilinear coordinates for multiple scattering[END_REF][START_REF] Alzubaidi | Formulation and accuracy of On-Surface Radiation Conditions for acoustic multiple scattering problems[END_REF][START_REF] Bruno | A new high-order high-frequency integral equation method for the solution of scattering problems. II: Multiple-scattering configurations[END_REF][START_REF] Chen | Scattering of sound from point sources by multiple circular cylinders using addition theorem and superposition technique[END_REF][START_REF] Ehrhardt | Wave Propagation in Periodic Media Analysis, Numerical Techniques and practical Applications[END_REF][START_REF] Ehrhardt | Numerical simulation of waves in periodic structures[END_REF][START_REF] Geuzaine | On the O(1) solution of multiple-scattering problems[END_REF][START_REF] Grote | Dirichlet-to-Neumann boundary conditions for multiple scattering problems[END_REF][START_REF] Kharlamov | Generalisation of the method of images for the calculation of inviscid potential flow past several arbitrarily moving parallel circular cylinders[END_REF][START_REF] Thierry | mu-diff: An open-source Matlab toolbox for computing multiple scattering problems by disks[END_REF][START_REF] Van Genechten | A Trefftz-based numerical modelling framework for Helmholtz problems with complex multiple-scatterer configurations[END_REF][START_REF] Yilmaz | An iterative procedure for the diffraction problem of water waves by multiple cylinders[END_REF] where methods are available to solve scattering by complex structures.

In the present paper, we consider the two-dimensional case where we have M sound-soft, sound-hard or homogeneous dielectric/penetrable circular cylinders. Even if the shape of the obstacles is simple, numerically solving this problem for many obstacles and for large frequencies remains nontrivial. In particular, getting efficient and accurate computational solutions if of utmost importance for important applications (electromagnetics, optics, nanophotonics for instance) where a huge number of diffraction objects are considered to model structured or disordered media (see e.g. [START_REF] Bidaud | Plasmon-based nanolenses assembled on a well-defined DNA template[END_REF][START_REF] Doyle | Modeling the permittivity of two-phase media containing monodisperse spheres: Effects of microstructure and multiple scattering[END_REF][START_REF] Doyle | Simulation of elastic wave scattering in cells and tissues at the microscopic level[END_REF][START_REF] Ferrand | Direct imaging of photonic nanojets[END_REF][START_REF] Hu | Compact wavelength demultiplexer via photonic crystal multimode resonators[END_REF][START_REF] Joannopoulos | Photonic Crystals: Molding the Flow of Light[END_REF][START_REF] Kharlamov | Generalisation of the method of images for the calculation of inviscid potential flow past several arbitrarily moving parallel circular cylinders[END_REF][START_REF] Natarov | Periodicity-induced effects in the scattering and absorption of light by infinite and finite gratings of circular silver nanowires[END_REF][START_REF] Savo | Walk dimension for light in complex disordered media[END_REF][START_REF] Tsang | Scattering of Electromagnetic Waves, Numerical Simulation[END_REF][START_REF] Wiersma | Disordered photonics[END_REF]).

To solve this class of problem, we first consider in Section 2 how to formulate the problem after writing the boundary-value problem. The formulation is developed following the Mie series expansion theory since the obstacles are circular cylinders. We then use the Dirichlet or Neumann boundary conditions to write the set of equations to solve. These equations are next suitably truncated and we explain quickly how to numerically solve them by spectral methods. We report some numerical examples (scattered/total fields and Radar Cross Sections) by using the Matlab toolbox µ-diff dedicated to solve such problems. Then we extend the methodology in Section 3 by using integral equation techniques in the case of dielectric circular cylinders. A numerical example is provided to illustrate the method. We end by a short conclusion in Section 4.

2 Multiple scattering by a cluster of M circular cylinders: Dirichlet and Neumann boundary conditions

Formulation

Let us consider a homogeneous medium for the whole space R 2 in which M disjoint scatterers Ω - 1 , . . . , Ω - M are included. We suppose that each of the scatterers Ω - p , p = 1, . . . , M , is a bounded domain in R 2 with boundary Γ p := ∂Ω - p . We define Ω -= ∪ M p=1 Ω - p as the global domain built from these separated obstacles. To start, we propose to analyze the case of the scattering problem of an incident plane wave u inc (x) = e ikβ•x (with x = (x, y) ∈ R 2 ), with direction β, by Ω -(assuming that the time dependence has the form e -iωt and for a real-valued wavenumber k). Concretely, this means that we wish to solve the following boundary-value problem: compute the scattered field u solution to (P)

       ∆u + k 2 u = 0, in Ω + := R 2 \ Ω -, γ(u) = -γ(u inc ), on Γ := ∂Ω -, lim ||x||→+∞ ||x|| 1/2 (∇u • x ||x|| -iku) = 0.
The operator ∆ = ∂ 2 x + ∂ 2 y is the Laplace operator and (∆ + k 2 ) is the Helmholtz operator. The gradient operator is ∇ and ||x|| = √ x • x, where x•y is the scalar product of two vectors x and y of R 2 . In practice, for sound-soft obstacles, the boundary operator γ is given by the zeroth-order trace operator, i.e. the boundary condition writes: u = -u inc on Γ. In the case of sound-hard obstacles, then the Neumann trace is involved leading to ∂ n u = -∂ n u inc on Γ, where n is the unit normal vector outwardly directed to the bounded scatterers Ω - p . The last equation is the well-known Sommerfeld's radiation condition at infinity which is needed to prove the existence of the solution [START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF][START_REF] Nédélec | Acoustic and Electromagnetic Equations[END_REF].

One possible and physically direct way to try solving the multiple scattering problem (P) is to interpret it as a collection of M single scattering problems that are next coupled through multiple scattering phenomena. This means that we can rewrite equivalently the problem by introducing the scattered fields u 1 , . . . , u M , considering that each field u p is given as the wave reflected only by the obstacle p illuminated by the combination of the incident wave and the M -1 scattered waves u q , for q = 1, . . . , M , with q = p.

In fact it can be proved [START_REF] Balabane | Boundary decomposition for Helmholtz and Maxwell equations. I. Disjoint subscatterers[END_REF] that the following result holds. We consider u as the solution to the multiple scattering problem (P). Then, we can show that the M coupled scattering problems, for p = 1, . . . , M , given by

(P p )              ∆u p + k 2 u p = 0, in R 2 \ Ω - p , γ(u p ) = -γ(u inc + M q=1,q =p u q ) on Γ p , lim ||x||→+∞ ||x|| 1/2 (∇u • x ||x|| -iku) = 0,
have a unique solution given by (u 1 , . . . , u p ). In addition, we also have the following decomposition of the wave field

u = M p=1 u p . (1)
We now consider the specific case where the scatterers Ω - p are circular cylinders of radius a In the sequel, we will also need the local polar coordinates of P in the orthonormal system of coordinates associated to the p-th scatterer, i.e. r p = --→ O p P, r p = |r p |, θ p = Angle(Ox, r p ). We define, for all m ∈ Z, the following cylindrical wave functions that are particular solutions to the Helmholtz equations for r > 0

ψ m (r) = H (1) m (kr)e imθ , ψ m (r) = J m (kr)e imθ . (2)
In the above relations, J n is known as the n-th order Bessel function while H

n is the n-th order first-kind

Hankel function [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF]. For all m ∈ Z, it can be shown that ψ m satisfies the outgoing Sommerfeld radiation condition. We also need, for all m ∈ Z, the local cylindrical wave functions related to the p-th scatterer,

for p = 1, . . . , M , ψ p m (r) = ψ m (r p ) = H (1) m (kr p )e imθp , ψ p m (r) = ψ m (r p ) = J m (kr p )e imθp , ∀ m ∈ Z. (3)
Because u p is an outgoing solution to a single scattering problem outside a circular cylinder, we have the decomposition

u p (r) = m∈Z c p m ψ p m (r), ∀ p = 1, . . . , M, ∀r p > a p . (4)
The complex coefficients (c p m ) m∈Z are obtained by simply using the boundary condition on Γ p

γ(u p ) = -γ(u inc ) - M q=1,q =p γ(u q ). (5)
To get a more explicit expression of this relation, we need to write the fields u inc and u q , for q = p, in the local system of coordinates of the scatterer p. If β = (cos β, sin β), then we have [START_REF] Martin | Multiple Scattering, Interaction of Time-Harmonic Waves with N Obstacles, Encyplopedia of Mathematics and its Applications[END_REF]

u inc (r) = m∈Z d p m ψ p m (r), (6)
where

d p m = e ikβ•bp e im( π 2 -β) .
In addition, the separation theorem [START_REF] Martin | Multiple Scattering, Interaction of Time-Harmonic Waves with N Obstacles, Encyplopedia of Mathematics and its Applications[END_REF] implies that, for 1 ≤ p, q ≤ M , with p = q, we have

ψ q m (r) =        n∈Z S mn (b pq ) ψ p n (r) for r p < b pq , n∈Z S mn (b pq )ψ p n (r) for r p > b pq , ∀ m ∈ Z, (7)
setting S mn (b pq ) = ψ m-n (b pq ) and S mn (b pq ) = ψ m-n (b pq ). The infinite dimensional matrices S p,q = (S mn (b pq )) m,n∈Z and S p,q = ( S mn (b pq )) m,n∈Z are the so-called separation matrices. Based on (1), ( 4), ( 5) and [START_REF] Bidaud | Plasmon-based nanolenses assembled on a well-defined DNA template[END_REF], we obtain the infinite linear systems

C p + D p M q=1,q =p (S p,q ) T C q = B p ∀ p = 1, . . . , M, (8)
where C p = (c p n ) n∈Z is the infinite vector containing the coefficients of the cylindrical decomposition (4) of u p , (S p,q ) T denotes the transpose of the separation matrix S p,q between the obstacles B p and B q defined by S p,q = (S p,q mn ) m∈Z,n∈Z and S p,q mn = ψ m-n (b pq ), D p = (D p mn ) mn∈Z is the diagonal infinite matrix, with diagonal terms

D p m,m =          J n (ka p ) H (1) n (ka p )
for sound-soft obstacles,

J n (ka p ) H (1) n (ka p )
for sound-hard obstacles,

B p = -D p d p , where d p = (d p m ) m∈Z is the infinite dimensional vector with coefficients (6) of the cylindrical incident field decomposition. Another way to write [START_REF] Bruno | A new high-order high-frequency integral equation method for the solution of scattering problems. II: Multiple-scattering configurations[END_REF] is

AC = B ( 10
)
where, if I is the identity operator on 2 (C),

A =       I D 1 S 1,2 T . . . D 1 S 1,M T D 2 S 2,1 T I . . . D 2 S 2,M T . . . . . . D M S M,1 T D M S M,2 T . . . I       , C =      C 1 C 2 . . . C M      , B =      B 1 B 2 . . . B M      . ( 11
)
Let us remark that the wave field u can be computed locally to Ω p (for r p < min 1≤q≤M q =p b pq ) by

u(r) = m∈Z c p m ψ p m (r) + m∈Z ( M q=1,q =p n∈Z S nm (b pq )c q n ) ψ p m (r), (12)
while the scattering amplitude a(θ) such that

u(r) = e ikr √ r a(θ) + O 1 r , for r → +∞, is a(θ) = e -iπ/4 2 πk p=M p=1 e -ibpk cos(θ-αp) ( n∈Z e in(θ-π 2 ) c p n ). (13)
In addition, expressions for the traces and normal derivative traces can be obtained with respect to the boundary conditions.

Computational aspects

For a practical implementation, the linear (8) must be truncated. Let us define the number of finite modes for the truncation of the solution u p by 2N p + 1, where N p ∈ N, C p = (c p n ) n=-Np,...,Np as the finite vector with approximations of the first 2N p + 1 coefficients c p n of the cylindrical decomposition (6) of u p , S p,q is the (2N p + 1) × (2N q + 1) finite dimensional separation matrix involving only the interactions between the first modes of Ω - p and Ω - q such that S p,q = (S p,q mn ) -Np≤m≤Np,-Nq≤n≤Nq , S p,q mn = ψ m-n (b pq ),

AC = B, (14)
with

A ∈ C N,N (of size N = M p=1 (2N p + 1)) ( I p is the identity matrix of C 2Np+1 ) A =       I 1 D 1 S 1,2 T . . . D 1 S 1,M T D 2 S 2,1 T I 2 . . . D 2 S 2,M T . . . . . . D M S M,1 T D M S M,2 T . . . I M       , C =      C 1 C 2 . . . C M      , B =      B 1 B 2 . . . B M      . ( 15
)
In practice, a suitable choice of the truncation parameter is given by (see e.g. [START_REF] Thierry | mu-diff: An open-source Matlab toolbox for computing multiple scattering problems by disks[END_REF])

N p = ka p + 1 2 √ 2 ln(2 √ 2πka p ε -1 ) 2 3 (ka p ) 1/3 + 1 , (16)
setting [x] as the integer part of a real number x, and ε is the error bound on the Fourier coefficients.

Finally, the numerical solution to the linear system ( 14) is obtained by a direct solver for dense linear systems or a preconditioned Krylov subspace iterative solver [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF] combined with fast Toeplitz acceleration techniques [START_REF] Chen | Matrix Preconditioning Techniques and Applications[END_REF]. In this case, a low memory storage of the system can be considered. The choice mainly depends on the physical configuration [START_REF] Thierry | mu-diff: An open-source Matlab toolbox for computing multiple scattering problems by disks[END_REF].

Numerical examples

To illustrate the method, we present two numerical examples based on the Matlab toolbox µ-diff [START_REF] Thierry | mu-diff: An open-source Matlab toolbox for computing multiple scattering problems by disks[END_REF].

Similar experiments can be reproduced based on available scripts that can be run on any computer. We consider the square box ] -8; 8[ 2 in which we randomly distribute M = 200 circular scatterers with radii such that 0.1 ≤ a p ≤ 0.2, for p = 1, ...M . The incident field is a plane wave with wavenumber k = 10 and angle of incidence β = 180 degrees. Let us remain that the circles are not sticky and that a minimal distance of 0.01 is fixed here between their boundary to avoid such a situation. We report on Figure 1 the results. On Figure 1a, we represent the scattering configuration for the M = 200 disks. On plots 1b and 1c, we draw the amplitude of the scattered and total wave fields, respectively. We can observe the strong interaction between the scatterers (with of course no penetration of the field inside the obstacles).

Finally, we give the bistatic Radar Cross Section (RCS) in the direction θ on Figure 1d where a is given by ( 13). Similarly, with the same parameters (but for a different location of the scatterers), we report on Figures 2a-2d the case of a Neumann problem which corresponds to the soundhard case. We observe that the physical behavior of the wave field is strongly dependent on the kind of boundary condition that is involved. In particular, the Dirichlet boundary conditions reflects the field mainly in the forward direction while the field penetrates inside the cluster for the Neumann case, generating some important backscattering propagation. 3 Multiple scattering by a collection of penetrable cylinders

Physical problem

We consider now M regular bounded and non intersecting circular dielectric ElectroMagnetic (EM)

scatterers Ω - p , p = 1, ..., M . They are assumed to be distributed in R 2 , with boundary Γ p := ∂Ω - p . The global dielectric scatterer Ω -is then built as the collection of the M separate single dielectric obstacles, i.e.

Ω -= ∪ M p=1 Ω - p , with boundary Γ = ∪ M p=1 Γ p . The homogeneous isotropic exterior domain of propagation is Ω + = R 2 \ Ω -.
We assume that we have a time-harmonic incident plane wave u inc (x) = e ik + β•x illuminating Ω -, with incidence direction β = (cos β, sin β) and time dependence e -iωt , where ω is the wave pulsation. The exterior wave number is k + = ω √ ε 0 µ 0 while the interior wavenumbers are such that k - p = ω √ ε p µ p , defining (ε 0 , µ 0 ) as (respectively (ε p , µ p )) the electric permittivity and magnetic permeability in the vacuum (respectively in the obstacle Ω - p ). Let us introduce the wavenumber k defined as piecewise constant with value k + out of the single obstacles and k -inside. Then, the problem writes:

find the scattered field u + and the transmitted wave u -solution to the transmission problem

               ∆u -+ (k -) 2 u -= 0, in Ω -, ∆u + + (k + ) 2 u + = 0, in Ω + , u + -u -= -u inc , on Γ, ∂ n u + -χ∂ n u -= -∂ n u inc , on Γ, lim ||x||→+∞ ||x|| 1/2 ∇u + • x ||x|| -ik + u + = 0. (17)
The total physical field u tot writes: u tot = u + + u inc , outside and u tot = u -inside the obstacles. In the so-called Transverse-Magnetic mode (or Transverse-Electric mode), the unknown u := u ± stands for the z-component H z of the magnetic field H and the transmission condition has to be read with χ p = ε 0 ε -1 p , on the boundary Γ p , while for TE modes (or TM mode), then u is linked to the electric field E through u = E z and the transmission condition is given by χ p = µ 0 µ -1 p . We define ε| Γp = ε p , µ| Γp = µ p , and χ| Γp = χ p , accordingly to the polarization.

An integral equation formulation of the scattering problem (17)

For system [START_REF] Geuzaine | On the O(1) solution of multiple-scattering problems[END_REF], we will use the theory of integral equations to write the set of equations to numerically resolve [START_REF] Antoine | An introduction to operator preconditioning for the fast iterative integral equation solution of time-harmonic scattering problems[END_REF][START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF][START_REF] Nédélec | Acoustic and Electromagnetic Equations[END_REF]. This can be done in the case of general obstacles thanks to the theory of potential.

Let us remark that such an approach could also be used for Dirichlet, Neumann and even impedance boundary conditions.

Let G be the two-dimensional free-space Green's function defined by ∀x, y ∈ R 2 , x = y, G(x, y) = i 4 H

(1) 0 (k xy ).

Integral equations are essentially based upon the Helmholtz integral representation formula [START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF]. Indeed, if v is a solution to the Helmholtz equation in the unbounded connected domain Ω + and satisfies the Sommerfeld radiation condition, then we have the following relation

Γ -G(x, y)∂ n v(y) + ∂ ny G(x, y)v(y) dΓ(y) = v(x) if x ∈ Ω + , 0 otherwise. ( 18
)
If v -is now solution to the Helmholtz equation in the bounded domain Ω -, then we have

Γ -G(x, y)∂ n v -(y) + ∂ ny G(x, y)v -(y) dΓ(y) = 0 if x ∈ Ω + , -v -(x) otherwise. ( 19
)
The integrals on Γ has to be interpreted in the sense of the inner product in L 2 (Γ) Then we can express v and v -based on ( 18)- [START_REF] Gumerov | Multiple scattering from N spheres using multipole reexpansion[END_REF] as

f, g H -1/2 ,H 1/2 = Γ f gdΓ,
   v(x) = -L(∂ n v + | Γ )(x) -M(v + | Γ )(x), ∀x ∈ Ω + , v -(x) = L(∂ n v -| Γ )(x) + M(v -| Γ )(x), ∀x ∈ Ω -.
Furthermore, the single-and double-layer integral operators provide some outgoing solutions to the Helmholtz equation [START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF]. Indeed, for any densities ρ ∈ H -1/2 (Γ) and λ ∈ H 1/2 (Γ), Lρ and Mλ, and any linear combination of both these functions yields to outgoing solutions to the Helmholtz equation in R 2 \Γ for some boundary conditions. Let us now recall the jump relations [START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF]. For any x in Γ, the trace and normal derivative traces of the operator L are given by (the signs means that z tends towards

x from the exterior or the interior of Γ)

lim z∈Ω ± →x Lρ(z) = Lρ(x), lim z∈Ω ± →x ∂ nz Lρ(z) = ∓ 1 2 I + N ρ(x), (20)
where I is the identity operator, for x ∈ Γ, and where

Lρ(x) = Γ G(x, y)ρ(y)dΓ(y), N ρ(x) = Γ ∂ nx G(x, y)ρ(y)dΓ(y).
In the paper, we denote the boundary integral operators by a roman letter (e.g. L). To solve [START_REF] Geuzaine | On the O(1) solution of multiple-scattering problems[END_REF] for general smooth scatterers, we can consider the single-layer representation of u + and u -

u + = L + ρ + and u -= L -ρ -, (21)
where L + (resp. L -) is the single-layer operator associated with k + (respectively k -). The unknown (ρ + , ρ -) is next solution to the integral equation

L + -L - I 2 -N + χ( I 2 + N -) ρ + ρ - = -u inc ∂ n u inc , (22)
that we rewrite

Aρ := A 11 A 12 A 21 A 22 ρ + ρ - = b := b 1 b 2 . ( 23
)
The plus and minus superscripts in L and N refers to as the wave numbers k + and k -. Since Ω -= ∪ M p=1 Ω - p is multiply connected, all the integral operators can be written as blocks structured systems and ( 23)

reads      A 1,1 A 1,2 . . . A 1,M A 2,1 A 2,2 . . . A 2,M . . . . . . . . . . . . A M,1 A M,2 . . . A M,M           ρ 1 ρ 2 . . . ρ M      =      b 1 b 2 . . . b M      , (24)
where A p,q ρ q = (A q ρ q )| Γp . After some manipulations and using a Schur complement then one gets ρ + (assuming that we are not on an interior resonance of one of the dielectric cylinders)

A S,+ ρ + := (L + + χ -1 L -( I 2 + N -) -1 ( I 2 -N + ))ρ + = b S,+ := -u inc + χ -1 L -( I 2 + N -) -1 ∂ n u inc , (25)
and ρ -is obtained by: ρ

-= A S,-ρ + + b S,-:= (L -) -1 (L + ρ + + u inc ).
In the case of circular cylinders, then a numerical approach as the previous spectral one in Section 2 allows to simplify the expression of the integral operators and to solve the linear system efficiently by a direct or iterative solver [START_REF] Thierry | mu-diff: An open-source Matlab toolbox for computing multiple scattering problems by disks[END_REF].

Numerical example

We now provide on Figure 3 an example of computation based on the use of µ-diff. We consider M = 30 dielectric scatterers (0.5 ≤ a p ≤ 1, ε p = 0.2ε 0 , µ p = µ 0 , p = 1, ..., M ) illuminated by a plane wave of incidence (-1, 0), with wave number k = 20. On plot 3a, we report the amplitude of the total field and see that the wave inside the scatterers propagates with a lower frequency, then producing a complex scattering pattern behind the cluster. In addition, we report the RCS on Figure 3b showing that there is an energy peak in the forward direction.

Conclusion

We explained in this paper how to efficiently solve scattering by a cluster of separated circular cylinders.

The methodology is relatively general and can handle many kinds of boundary conditions. In particular, using integral equations provides the access to a wide variety of formulations. We did not developed much the explanations according to the numerical schemes since they are available in the Matlab toolbox µ-diff. Because we are using the theory of integral equations, all the current developments can be directly applied to the scattering by spheres for three-dimensional problems, in particular for acoustic waves even if similar theories exist for the EM and elastic cases [START_REF] Antoine | An introduction to operator preconditioning for the fast iterative integral equation solution of time-harmonic scattering problems[END_REF]. However, the numerical schemes must be more deeply investigated since the structure of the matrices defining the linear system to resolve strongly depends on the spherical harmonics expansions [START_REF] Gumerov | Multiple scattering from N spheres using multipole reexpansion[END_REF].

  p and centered at O p = (x p , y p ) of a given orthonormal system of coordinates ( Let us introduce the following notations: for all p = 1, . . . , M : b p = --→ OO p , b p = |b p |, α p = Angle(Ox, b p ) and for all q = 1, . . . , M , with q = p: b pq = ---→ O q O p , b pq = |b pq |, α pq = Angle(Ox, b pq ). Let us consider a point P of the plane described by its cartesian coordinates (x, y), or alternatively by its polar coordinates: r = --→ OP, r = |r|, θ = Angle(Ox, r).

  and D p = (D p mn ) -Np≤m≤Np,-Nq≤n≤Nq is the diagonal finite matrix, with diagonal terms D p m,m = D p m,m , B p = -D p d p , where d p = (d p m ) -Np≤m≤Np contains the 2N p + 1 first coefficients of the cylindrical incident field decomposition (4). The M coupled finite dimensional systems writes

  which allows to analyze where the main lobs related to the diffraction problem arise. Its definition is given by RCS(θ) := 10 log 10 (|a(θ)| 2 ) (dB),

Figure 1 :

 1 Figure 1: Scattering by M = 200 sound-soft disks randomly distributed in the box ] -8; 8[ 2 . The incident plane is fixed by k = 10 and β=(-1,0).

Figure 2 :

 2 Figure 2: Scattering by M = 100 sound-hard disks randomly distributed in the box ]-8; 8[ 2 . The incident plane is fixed by k = 2π and β=(-1,0).

  because both u inc and Γ are regular. Let us define the volume single-and double-layer integral operators, respectively denoted by L and M, as, ∀x ∈ R 2 \Γ, L : ρ → Lρ(x) = Γ G(x, y)ρ(y) dΓ(y), M : λ → Mλ(x) = -Γ ∂ ny G(x, y)λ(y) dΓ(y).

Figure 3 :

 3 Figure 3: Scattering by M = 30 dielectric disks randomly distributed in the box ] -8; 8[ 2 . The incident plane is fixed by k = 20 and β=(-1,0).

Acknowledgments. The authors would like to thank the Deanship of Scientific Research at Umm Al-Qura University for supporting this work by grant code 18-SCI-1-01-0017.