Elsa Piollet
email: elsa.piollet@polymtl.ca

Guilhem Michon
email: guilhem.michon@isae.fr

Annie Ross
email: annie.ross@polymtl.ca

TRANSIENT RESPONSE OF A SANDWICH STRUCTURE DAMPED WITH A FIBROUS CORE MATERIAL

Keywords: Vibration damping, sandwich structures, nonlinear behavior, transient response

This paper investigates the effect of a fibrous core material on the transient response of a sandwich structure. The core material studied is made of entangled carbon fibers crosslinked with epoxy resin. In order to understand its complex behavior, the material is characterized in shear oscillation, and the obtained shear-stress loops are described using Dahl's dynamic hysteresis model. Then, the transient response of a single-degree-of-freedom system including this material is simulated and the amplitude dependency of the response is analyzed. Sandwich beams with the fibrous core material are tested experimentally with a impact hammer and compared with classical honeycomb and foam cored sandwich beams. The entangled cross-linked fibers are shown to provide high damping resulting in a fast return to the equilibrium position. Both theoretical and experimental studies showed nonlinear damping. In the amplitude range studied, the material is more interesting at high impact amplitude.

INTRODUCTION

Aerospace and automotive structures need to be as light as possible in order to limit fuel energy consumption, while maintaining high stiffness and strength properties. Sandwich structures, composed of two facesheets separated by a light core, are widely used in such applications as they provide high stiffness-to-weight ratios. However, the use of stiff lightweight structures can lead to noise and vibration amplification. This can be detrimental to user comfort and structural durability. One way to reduce noise and vibration is to investigate new core materials with high energy dissipation properties.

Studies on the static behavior of entangled fibers indicate that energy is dissipated by friction between the fibers when they are deformed [START_REF] Poquillon | Experimental data about mechanical behaviour during compression tests for various matted fibres[END_REF]. This has led to the development of a new sandwich core material in which entangled fibers are cross-linked with epoxy resin [START_REF] Mezeix | Mechanical behavior of entangled fibers and entangled crosslinked fibers during compression[END_REF]. Resin cross-links ensure cohesion and stiffness of the material while fibers can provide energy dissipation. The present study investigates the effect of this material on the damping properties of sandwich structures through a transient analysis.

In this paper, the material behavior is first measured and identified with a hysteresis model. Then, the transient response of a single-degree-of-freedom system is studied numerically. The consequences of the material behavior on this simple structural response is analyzed. Finally, the transient response of a sandwich beam with entangled core material is compared experimentally with the response of classical honeycomb and foam sandwich beams.

MATERIAL CHARACTERIZATION

The material is made following the fabrication process described by Mezeix [START_REF] Mezeix | Mechanical behavior of entangled fibers and entangled crosslinked fibers during compression[END_REF]. Different types of fibers can be used. Here, carbon fibers are chosen, since carbon-epoxy materials are commonly used in aerospace applications. The fibers are 7 µm in diameter, with a Young's modulus of 240 GPa and a density of 1770 kg/m 3 . They are first cut to a length of 31 mm, and are then entangled in a blower room under an air flux at a pressure of 5 bar.

Epoxy resin is then sprayed on the entangled fibers with The final density of the material is 180 kg/m 3 (150 kg/m 3 of carbon fibers and 30 kg/m 3 of epoxy resin). Figure 1 shows a scanning electron microscopy of the polymerized material. It can be seen that some contacts are bonded by resin while others are free to slip.

As the material is intended to be used as a core material in sandwich structures, shear is the most important mechanism to be analyzed [START_REF] Mead | The forced vibration of a three-layer, damped sandwich beam with arbitrary boundary conditions[END_REF]. Figure 2 shows the set-up for the shear test on a BOSE ElectroForce R 3330. Two samples of dimensions 20 mm × 40 mm × 60 mm are tested together to ensure shearing only. The samples are held by two steel brackets clamped on the bottom plate. A vertical displacement is imposed between the samples by an aluminum plate attached to the vertical motor. The resulting force is measured. From the displacement and force, the shear strain γ and the shear stress τ in the samples are deduced by assuming a constant shear through the thickness.

Measurements are made at 1 Hz, 20 Hz, 40 Hz, 60 Hz and 80 Hz, and the material behavior is found to be frequency independent in the tested frequency range. The following study is based on measurements made at 20 Hz, as they were found to be less sensitive to noise coming from the set-up.

The measured shear stress can be decomposed as a linear part of modulus G 1 and a hysteresis part defined as:

τ H = τ -G 1 γ (1)
The measured shear stress without the linear part, τ H , is represented against the measured shear strain γ on Fig. 3 (grey line). The physical interpretation of the observed hysteresis is the following: after each change in the strain direction, the strain rate is null, and all unbonded contacts between fibers are blocked in a "stick" configuration, which creates an initial stiffness; then, stiffness decreases as contacts begin to slip individually, until they are all in a "slip" configuration. As for contacts cross-linked with epoxy resin, they never enter the "slip" configuration, and they are responsible for the linear behavior G 1 γ which does not appear on the figure.

The shape of the hysteresis loop, combined with the frequency independent behavior, indicates that dissipation comes from the friction between the fibers.

The hysteresis loops are described using Dahl's dynamic hysteresis model [START_REF] Dahl | Math model of hysteresis in piezo-electric actuators for precision pointing systems[END_REF]:

dτ H dγ = σ τ C (2τ C -(τ H -τ m ) sgn( γ)) (2)
where γ is the strain rate, σ and τ C are constants of the model and τ m is the value of the stress at the last change in the strain direction.

Figure 3 shows the identified hysteresis loops for one set of samples (black dotted line). Two sets of two samples each are tested, leading to the following parameters for the model:

G 1 = 6.02 • 10 6 Pa σ = 1.49 • 10 6 Pa τ C = 1.39 • 10 3 Pa (3)

SINGLE-DEGREE-OF-FREEDOM SYSTEM

A theoretical single-degree-of-freedom system is now simulated in order to understand the consequences of the measured material properties on the transient response of a structure. The structure studied is composed of two samples in shear configuration, as tested previously, but separated by a rigid mass. Figure 4 shows the principle of the system. The mass is taken to be m = 1 kg, and is assumed to move only in the vertical direction. The parameters of the material are as given in Eq. (3). A base acceleration is applied at the bottom of the structure. The shock acceleration is a half sine of duration 1.5 ms, as shown on Fig. 5(a). Transient response is computed for seven amplitudes ranging from 1 g to 25 g.

Figure 5(b) shows the simulated transient response from 0 s to 0.2 s. The response to the 25 g shock shows three decay regimes. First, at high oscillation amplitudes, fast linear decay typical of friction is observed (regime 1). Then, the decay seems to become exponential, which is more representative of viscous damping. However, the slope decreases with the oscillation amplitude, from relatively fast (regime 2) to very slow decay (regime 3). Thus, while most of the oscillation is damped before 0.1 s, the remaining oscillation is dissipated slowly.

For lower shock amplitudes, the response enters directly the "exponential" decay regimes 2 and 3, as shown on Fig. 5(c) which spans 0.5 s. The resulting decay rate is slower for those low amplitude oscillations, and the remaining oscillations when entering regime 3 represent a large part of the initial oscillation amplitude.

This simulation shows that the material behavior leads to a transient response that depends on the shock amplitude. For this simple structure and in the amplitude range studied, the entangled material leads to higher damping for higher excitation levels. The following part presents the behavior of a real structure studied experimentally in order to see if the same phenomena are encountered.

TRANSIENT RESPONSE OF SANDWICH BEAMS

The transient responses of three sandwich beams are studied experimentally. For the purpose of comparing with standard materials, the core materials of the beams are the following, along with their density ρ and linear shear modulus G:

• entangled cross-linked carbon fibers (ρ = 180 kg/m 3 , G 1 = 6.02 MPa) • Nomex honeycomb (ρ = 48 kg/m 3 , G = 44.8 MPa)

Time (s) Displacement x (m) 0 0.05 0.1 0.15 0 .2 -2 -1 0 1 2 x 10 -4 (b)
• PMI foam (ρ = 52 kg/m 3 , G = 19 MPa).
The beam length is 253 mm and the width is 40 mm. The skins are made of seven plies of prepreg carbon-epoxy fabric used in the aeronautic industry, for a total thickness of 2 mm. The core thickness is 20 mm.

The set-up is shown on Fig. 6. The sandwich beams are tested in cantilever configuration. Excitation is applied with an impact hammer close to the clamped end, and the response is measured at the free end with a laser vibrometer.

Figure 7 shows the measured impacts and transient response for two impact amplitudes of 10 N and 100 N. The durations of the impacts were 1.7 ms and 1.4 ms respectively. The response shown was obtained as the average of five successive tests. These measurements show that the entangled cross-linked fibers lead to a large improvement of damping as compared to classical core materials. While the PMI foam dissipates more than the honeycomb, the entangled cross-linked fibers are far more dissipative. The difference increases with impact amplitude, and the material is thus more interesting for high impact amplitudes.

As expected from the previous simulation, the initial slope of the response is steep, particularly at 100 N. Then, the decay rate decreases, and at low oscillation amplitude the material provides little dissipation.

While this study demonstrates very good damping properties, it should be noted that the entangled cross-linked fibers tested have a higher density and a lower stiffness than the two other core materials.

CONCLUSIONS

The effect of a fibrous material on the transient response of a structure has been studied. Material study showed a hysteretic behavior indicating dissipation coming from friction between the fibers. Dahl's dynamic hysteresis model was used to model the measured hysteresis. The simulated transient response of a single-degree-of-freedom system showed three decay regimes depending on the amplitude. The first regime, typical of friction, corresponds to the higher decay rate, and is thus the most interesting. For the amplitude range studied, the entangled material is more dissipative for higher excitation amplitudes.

The transient response of a sandwich beam with entangled core material showed high added damping as compared to two classical core materials (Nomex honeycomb and PMI foam). Moreover, as predicted by the single-degree-of-freedom simulation, the added damping was increased at higher impact amplitude.

This study shows that the entangled-cross-linked material is a good candidate for vibration damping of sandwich structures. It should be noted however that its lower stiffness and higher density than classical core materials prevents a total replacement of classical materials. Future studies should be conducted to investigate possible applications and analyze the effect of the fiber density and fiber type on damping.

Figure 1 :

 1 Figure 1: Scanning electron microscopy of the entangled cross-linked fibers after polymerization.

Figure 2 :

 2 Figure 2: Experimental set-up for material characterization in shear.

Figure 3 :

 3 Figure 3: Shear stress-strain hysteresis loops obtained by removing the linear part of the shear stress: τ H = τ -G 1 γ with G 1 identified in Eq. (3). Loops measured at 20 Hz (grey line) and identified (black dotted line).

Figure 4 :

 4 Figure4: Principle of the studied single-degree-of-freedom system, including the entangled fiber samples in grey.

Figure 5 :

 5 Figure 5: Simulated transient response of the single-degreeof-freedom system of Fig. 4: (a) imposed acceleration, (b) response for all amplitudes during 0.2 s and (c) zoom on the low response amplitudes during 0.5 s. Numbers indicate decay regimes: linear decay (1), and fast (2) and slow (3) "exponential" decay.

Figure 6 :

 6 Figure 6: Experimental set-up for transient testing.

ForceFigure 7 :

 7 Figure 7: Transient response of sandwich beams with three different core materials to two impacts: (a) 10 N, 1.7 ms and (b) 100 N, 1.4 ms.

ACKNOWLEDGMENTS

This study was carried out during the first author's PhD work at Institut Supérieur de l'Aéronautique et de l'Espace (ISAE) in Toulouse, France. Midi Pyrénée region, Université de Toulouse and Ecole Doctorale Aéronautique et Astronautique are gratefully acknowledged for their financial support.


Adrien Clarenne
email: adrien.clarenne@univ-rennes1.fr

Rescaled weighted determinantal random balls

published or not. The documents may come

Rescaled weighted determinantal random balls

Introduction

In this work, we give a generalization of the existing results concerning the asymptotics study of random balls model. The first results are obtained in 2007 by Kaj, Leskela, Norros and Schmidt in [START_REF] Kaj | Scaling limits for random fields with long-range dependence[END_REF]. In their model, the balls are generated by an homogeneous Poisson point process on R d × R + (see [START_REF] Daley | Introduction to point processes[END_REF] for a general reference on point processes). In 2009, Breton and Dombry generalize this model adding in [START_REF] Breton | Rescaled weighted random balls models and stable self-similar random fields[END_REF] a mark m on the balls of the previous model and they obtain limit theorems on the so-called rescaled weighted random balls model. In 2010, Biermé, Estrade and Kaj obtain in [START_REF] Biermé | Self-similar random fields and rescaled random balls models[END_REF] results performing for the first time a zoom-in scaling. In 2014, Gobard in his paper [START_REF] Gobard | Random balls model with dependence[END_REF] extends the results of [START_REF] Breton | Rescaled weighted random balls models and stable self-similar random fields[END_REF] considering inhomogeneous weighted random balls, and adding a dependence between the centers and the radii. The next step is to consider repulsion between the balls. In [3], Breton, Clarenne and Gobard give results on determinantal random balls model, but no weight are considered in their model. In this note, we consider weighted random balls generated by a non-stationary determinantal point process. To that purpose, we give an extension of the Laplace transform of determinantal processes allowing to compute Laplace transform with not necessarily compactly supported function, but with instead a condition of integrability. The main contributions of this note thus are a simplication of the proof of [3] and the introduction of weights in the non-stationary determinantal random balls model. This paper is organized as follows. In Section 1 and Section 2, we give a description of the model and state our main results under the three different regimes. In Section 3, we give the Laplace transform of a determinantal point process for not compactly supported test functions and prove our results. Finally some technical results are gathered in the Appendix.

Model

We consider a model of random balls in R d constructed in the following way. The centers of the balls are generated by a determinantal point process (DPP) φ on R d characterized by its kernel K with respect to the Lebesgue measure. The motivation for considering such processes is that it introduces repulsion between the centers in agreement with various real model of balls. We assume that the map K given for all f ∈ L 2 (R d , dx) and x ∈ R d by

Kf (x) = R d K(x, y)f (y) dy (1)
satisfies the following hypothesis:

Hypothesis 1 The map K given in (1) is a bounded symmetric integral operator K from L 2 (R d , dx) into L 2 (R d , dx), with a continuous kernel K with spectrum included in [0, 1[. Moreover, K is locally trace-class, i.e. for all compact Λ ⊂ E, the restriction K Λ of K on L 2 (Λ, λ) is of trace-class.
Moreover, we also assume

x -→ K(x, x) ∈ L ∞ (R d ). (2)
These assumptions imply that K(x, x) ≥ 0. At each center x ∈ R d , we attach two positive marks r and m independently. The first mark is interpreted as the radius and the second mark is the weight of the ball B(x, r).

The radii (resp. the weight) are independently and identically distributed according to F (resp. according to G), assumed to admit a probability density f (resp. a probability density g). We have a new point process Φ on

R d × R + × R + and according to Proposition A.7 in [3], Φ is a DPP on R d × R + × R + with kernel K (x, r, m), (y, s, m ′ ) = g(m) f (r)K(x, y) f (s) g(m ′ ),
with respect to the Lebesgue measure. Moreover, we suppose that the probability measure G belongs to the normal domain of attraction of the α-stable distribution S α (σ, b, τ

) with α ∈ (1, 2]. Because α > 1, we can note that R + mG(dm) = R + mg(m)dm < +∞. (3)
In the sequel, we shall use the notation Φ both for the marked DPP (i.e. the random locally finite collection of points (X i , R i , M i )) and for the associated random measure (X,R,M )∈Φ δ (X,R,M ) . We consider the contribution of the model in any suitable measure µ on R d given by the following measure-indexed random field:

M(µ) = R d ×R + ×R + mµ B(x, r) Φ(dx, dr, dm). (4)
However, in order to ensure that M(µ) in ( 4) is well defined, we restrain to measures µ with finite total variation (see below Proposition 1.1). In the sequel, Z(R d ) stands the set of signed (Borelian) measures µ on R d with finite total variation µ var (R d ) < +∞. Moreover as in [START_REF] Kaj | Scaling limits for random fields with long-range dependence[END_REF], we assume the following assumption on the radius behaviour, for d < β < 2d,

f (r) ∼ r→+∞ C β r β+1 , r β+1 f (r) ≤ C 0 . (5)
Since β > d, condition [START_REF] Daley | Introduction to point processes[END_REF] implies that the mean volume of the random ball is finite:

v d +∞ 0 r d f (r) dr < +∞, (6)
where v d is the Lebesgue measure of the unit ball of R d . On the contrary, β < 2d implies that F does not admit a moment of order 2d and the volume of the balls has an infinite variance. The asymptotics condition in ( 5) is of constant use in the following. Proof: The proof follows the same lines as that of Proposition 1.1 in [3], replacing K(0) by K(x, x) and controlling it thanks to Hypothesis (2).

Asymptotics and main results

The zooming-out procedure acts accordingly both on the centers and on the radii. First, a scaling S ρ : r → ρr of rate ρ ∈ (0, 1] changes balls B(x, r) into B(x, ρr); this scaling changes the distribution F of the radius into

F ρ = F • S -1 ρ .
Second, the intensity of the centers is simultaneously adapted; to do this, we introduce actually a family of new kernels K ρ , ρ ∈]0, 1], that we shall refer to as scaled kernels, and we denote by φ ρ the DPP with kernel K ρ (with respect to the Lebesgue measure).The zoom-out procedure consists now in introducing the family of DPPs φ ρ , ρ ∈]0, 1], with kernels K ρ with respect to the Lebesgue measure satisfying

K ρ (x, x) ∼ ρ→0 λ(ρ)K(x, x), (7)
with lim ρ→0 λ(ρ) = +∞. We also suppose sup

x∈R d K ρ (x, x) ≤ λ(ρ) sup x∈R d K(x, x), (8)
and observe that with ( 2) and ( 8), Proposition A.6 in [3] gives the following uniform bound sup

x∈R d R d K ρ (x, y) 2 dy = ρ→0 O λ(ρ) . (9)
The zoom-out procedure consists in considering a new marked DPP Φ ρ on R d × R + × R + with kernel:

K ρ (x, r, m), (y, s, m ′ ) = g(m) f (r/ρ) ρ K ρ (x, y) f (s/ρ) ρ g(m ′ ),
with respect to the Lebesgue measure. The so-called scaled version of M(µ) is then the field

M ρ (µ) = R d ×R + ×R + mµ B(x, r) Φ ρ (dx, dr, dm).
In the sequel, we are interested in the fluctuations of M ρ (µ) with respect to its expectation

E M ρ (µ) = R d ×R + ×R + mµ B(x, r) K ρ (x, x) f (r/ρ) ρ g(m)dxdrdm
and we introduce

M ρ (µ) = M ρ (µ) -E M ρ (µ) = R d ×R + ×R + mµ B(x, r) Φ ρ (dx, dr, dm), (10)
where Φ ρ stands for the compensated random measure associated to Φ ρ .

We introduce a subspace M α,β ⊂ Z on which we will investigate the convergence of the random field M ρ (µ). The next definition comes from [START_REF] Breton | Rescaled weighted random balls models and stable self-similar random fields[END_REF]. Definition 2.1 For 1 < α ≤ 2 and β > 0, we denote by M α,β the subset of measures µ ∈ Z(R d ) satisfying for some finite constant C µ and some 0 < p < β < q :

R d |µ(B(x, r))| α dx ≤ C µ (r p ∧ r q )
where a ∧ b = min(a, b).

We denote by M + α,β the space of positive measures µ ∈ M α,β . Now, we can state the main result of this note. The proof consists in a combination of the arguments of [3] and [START_REF] Gobard | Random balls model with dependence[END_REF]. It is given in Section 3 where for some required technical points, it is referred to [3] and [START_REF] Breton | Rescaled weighted random balls models and stable self-similar random fields[END_REF].

Theorem 2.2 Assume (5) and φ ρ is a DPP with kernel satisfying (2), ( 7), ( 8) and Hypothesis 1 for its associated operator K ρ in [START_REF] Hough | Zeros of Gaussian Analytic Functions and Determinantal Point Processes[END_REF].

(i) Large-balls scaling: Assume λ(ρ)ρ β → +∞ and set n(ρ) = λ(ρ)ρ β 1/α . Then, M ρ (•)/n(ρ) converges in the fdd sense on M + α,β to W α (•)
where

W α (µ) = R d ×R + µ B(x, r) M α (dx, dr)
is a stable integral with respect to the α-stable random measure M α with control measure σ α K(x, x)C β r -β-1 dxdr and constant skewness function b given in the domain of attraction of G.

(ii) Intermediate scaling: Assume λ(ρ)ρ β → a d-β ∈]0, +∞[ and set n(ρ) = 1. Then, M ρ (•)/n(ρ) converges in the fdd sense on M + α,β to P • D a where P (µ) = R d ×R + ×R + mµ B(x, r) Π(dx, dr, dm)
with Π a (compensated) PPP with compensator measure K(x, x)C β r -β-1 dxdrG(dm) and D a is the dilatation defined by (D a µ)(B) = µ(a -1 B).

(iii) Small-balls scaling: Suppose λ(ρ)ρ β → 0 when ρ → 0 for d < β < αd and set

n(ρ) = (λ(ρ)ρ β ) 1/γ with γ = β/d ∈]1, α[. Then, the field n(ρ) -1 M ρ (•) converges when ρ → 0 in the finite-dimensional distributions sense in L 1 (R d ) ∩ L 2 (R d ) ∩ {µ ≥ 0} to Z γ (•) where Z γ (µ) = R d φ(x) M γ (dx) for µ(dx) = φ(x)dx with φ ∈ L 1 (R d ) ∩ L 2 (R d ), φ ≥ 0,
is a stable integral with respect to the γ-stable random measure M γ with control measure σ γ K(x, x)dx where

σ γ γ = C β v γ d d +∞ 0 1 -cos(r) r 1+γ dr +∞ 0 m γ G(dm),
and constant unit skewness.

Here, and in the sequel, we follow the notations of the standard reference [START_REF] Samorodnitsky | Stable Non-Gaussian Random Processes[END_REF] for stable random variables and integrals.

Proof

To investigate the behaviour of M ρ (µ) in the determinantal case, we use the Laplace transform of determinantal measures. An explicit expression is well known when the test functions are compactly supported, see Theorem A.4 in [3]. However, in our situation, the test functions (x, r, m) -→ mµ(B(x, r)) are not compactly supported on

R d × R + × R + for µ ∈ M + α,β .
In order to overpass this issue we use Proposition 3.1 below for the Laplace transform of determinantal measures with non-compactly supported test functions, but with a condition of integration with respect to the kernel of the determinantal process (see ( 11)). In addition to generalizing the model studied in [3] by adding a weight, the following proposition has the further consequence of simplifying the proofs of the results in [3], since there is no more need to study the truncated model and obtain uniform convergence to exchange the limit in R, the truncation parameter and the limit in ρ, the scaling parameter.

Proposition 3.1 Let Φ a determinantal point process on a locally compact Polish space E with a continuous kernel K such that the associated operator K satisfies Hypothesis 1. We also suppose that the kernel

K 1 -e -h ∈ L 2 (E × E). Let h be a nonnegative function such that E 1 -e -h(x) K(x, x)dx < +∞. ( 11
)
Then K 1 -e -h is a trace-class operator with

Tr K 1 -e -h = E 1 -e -h(x) K(x, x)dx
and we have

E exp - E h(x)Φ(dx) = exp - ∞ n=1 1 n Tr K 1 -e -h n , ( 12
)
where K 1 -e -h is the operator with kernel y) .

K 1 -e -h (x, y) = 1 -e -h(x) K(x, y) 1 -e -h(
Proof: Expression (12) is known to be true when h has a compact support (see Theorem A.4 and equation (37) in [3]), but it is not the case here. Let (h p ) p∈N a non-decreasing sequel of positive functions with compact support defined by

h p (x) = h(x)1 B(0,p) (x).
Thanks to Theorem A.4 in [3], we have for all p ∈ N:

E exp - E h p (x)Φ(dx) = exp - ∞ n=1 1 n Tr K 1 -e -hp n . ( 13
)
First step: To begin, we prove that

lim p→+∞ E exp - E h p (x)Φ(dx) = E exp - E h(x)Φ(dx) .
If we denote by M p (E) the space of all point measures defined on E, we have

E exp - E h p (x)Φ(dx) = Mp(E) exp - E h p (x)m(dx) P Φ (dm).
Because p -→ h p is increasing, by monotone convergence we have

lim p→+∞ E h p (x)m(dx) = E h(x)m(dx).
To finish we apply the dominated convergence theorem:

• lim p→+∞ exp - E h p (x)m(dx) = exp - E h(x)m(dx) ,
• exp -E h p (x)m(dx) ≤ 1 because h p ≥ 0, and 1 is integrable on M p (E) with respect to P Φ .

Second step: We prove that K 1 -e -h is a trace-class operator.

Let (e n ) n∈N be an orthonormal basis of L 2 (E).

K 1 -e -h is trace-class if +∞ n=0 K 1 -e -h (e n
), e n < +∞ and in this case we have

Tr K 1 -e -h = +∞ n=0 K 1 -e -h (e n ), e n .
We begin by proving that for all n ∈ N,

lim p→+∞ K 1 -e -hp (e n ), e n = K 1 -e -h (e n ), e n . K 1 -e -hp (e n ), e n = E×E
1 -e -hp(x) K(x, y) 1 -e -hp(y) e n (x)e n (y)dxdy

We apply the dominated convergence theorem:

• Computation of the limit when p → +∞: lim p→+∞ 1 -e -hp(x) K(x, y) 1 -e -hp(y) e n (x)e n (y) = 1 -e -h(x) K(x, y) 1 -e -h(y) e n (x)e n (y)

• Domination:

1 -e -hp(x) y) 1 -e -hp(y) e n (x)e n (y)

≤ 1 -e -h(x) |K(x, y)| 1 -e -h(y) |e n (x)| |e n (y)| which is integrable on E × E because: E×E 1 -e -h(x) |K(x, y)| 1 -e -h(y) |e n (x)| |e n (y)| dxdy ≤ E×E (1 -e -h(x) ) |K(x, y)| 2 (1 -e -h(y) )dxdy × E×E |e n (x)| 2 |e n (y)| 2 dxdy ≤ E (1 -e -h(x) )K(x, x)dx < +∞ because |K(x, y)| 2 ≤ K(x, x)K(y, y) and E |e n (x)| 2 dx = 1.
Thus we have:

+∞ n=0 K 1 -e -h (e n ), e n = +∞ n=0 lim inf p→+∞ K 1 -e -hp (e n ), e n ≤ (a) lim inf p→+∞ +∞ n=0 K 1 -e -hp (e n ), e n = lim inf p→+∞ Tr K 1 -e -hp
where (a) is allowed by the Fatou lemma because K 1 -e -hp (e n ), e n ≥ 0. We have

Tr K 1 -e -hp = E 1 -e -hp(x) K(x, x)dx.
We apply the dominated convergence theorem:

1. lim p→+∞ 1 -e -hp(x) K(x, x) = 1 -e -h(x) K(x, x),
2. For all p ≥ 0 we have:

1 -e -hp(x) K(x, x) ≤ 1 -e -h(x) K(x, x)
which is integrable on E by hypothesis (11).

Thus lim p→+∞ E

1 -e -hp(x) K(x, x)dx

= E 1 -e -h(x) K(x, x)dx. +∞ n=0 K 1 -e -h (e n ), e n ≤ E
1 -e -h(x) K(x, x)dx < +∞ which proves that K 1 -e -h is trace-class.

Third step: We prove that Tr K 1 -e -h = E

1 -e -h(x) K(x, x)dx.

Because K 1 -e -h is trace-class and K 1 -e -h ∈ L 2 (E×E), Lemma 4.2.2 in [START_REF] Hough | Zeros of Gaussian Analytic Functions and Determinantal Point Processes[END_REF] ensures that

1 -e -h(x) K(x, y) 1 -e -h(y) = +∞ i=1 λ i ϕ i (x)ϕ i (y) ( 14
)
where (ϕ i ) is an orthonormal basis of L 2 (E). Thus, thanks to (14) we have:

K 1 -e -h (e n ), e n = E×E 1 -e -h(x) K(x, y) 1 -e -h(y) e n (x)e n (y)dxdy = E×E +∞ i=1 λ i ϕ i (x)ϕ i (y)e n (x)e n (y)dxdy (15)
We can invert the sum and the integral if

+∞ i=1 E×E |λ i ϕ i (x)ϕ i (y)e n (x)e n (y)| dxdy < +∞. +∞ i=1 E×E |λ i ϕ i (x)ϕ i (y)e n (x)e n (y)| dxdy = +∞ i=1 E×E λ i |ϕ i (x)| |ϕ i (y)| |e n (x)| |e n (y)| dxdy = +∞ i=1 λ i E |ϕ i (x)| |e n (x)| dx 2 ≤ (1) +∞ i=1 λ i E |ϕ i (x)| 2 dx E |e n (x)| 2 dx = (2) +∞ i=1 λ i E |ϕ i (x)| 2 dx = (3) E +∞ i=1 λ i |ϕ i (x)| 2 dx = (4) E (1 -e -h(x) )K(x, x)dx < +∞.
The inequality ( 1) is the Cauchy-Schwarz inequality, the equality ( 2) is because

E |e n (x)| 2 dx = e n 2 2 = 1,
the inversion of the sum and integral in (3) is allowed because all the terms are nonnegative and the equality (4) stands because K is continuous. So we can invert the integral and the sum in (15) to obtain

K 1 -e -h (e n ), e n = +∞ i=1 E×E λ i ϕ i (x)ϕ i (y)e n (x)e n (y)dxdy = +∞ i=1 λ i E ϕ i (x)e n (x)dx 2 = +∞ i=1 λ i ϕ i , e n 2
and the result is proved because

+∞ n=0 K 1 -e -h (e n ), e n = +∞ n=0 +∞ i=1 λ i ϕ i , e n 2 = (a) +∞ i=1 λ i +∞ n=0 ϕ i , e n 2 = +∞ i=1 λ i ϕ i 2 2 = +∞ i=1 λ i E |ϕ i (x)| 2 dx = (a) E +∞ i=1 λ i |ϕ i (x)| 2 dx = E (1 -e -h(x) )K(x, x)dx < +∞,
where the two equalities (a) are allowed because the terms are non-negative. Then, K 1 -e -h is a trace-class operator with

Tr K 1 -e -h = E 1 -e -h(x) K(x, x)dx.
The computations above gives in particular

lim p→+∞ Tr K 1 -e -hp = Tr K 1 -e -h .
Fourth step: Now, we have to take the limit in the right term of equation ( 13).

we prove that we can exchange the limit and the infinite sum. To do that, we show that the sum normally converges. For all n, p ≥ 0:

Tr K 1 -e -hp n ≤ K 1 -e -hp n-1 Tr K 1 -e -hp ≤ K 1 -e -h n-1 Tr K 1 -e -h ,
where the second inequality stands thanks to Lemma 4.1 in Appendix 4 and because K(x, x) ≥ 0,

Tr K 1 -e -hp = E 1 -e -hp(x) K(x, x)dx ≤ E 1 -e -h(x) K(x, x)dx = Tr K 1 -e -h , (16)
which is finite as proved above. Let now prove that K 1 -e -h < 1.

K 1 -e -h = sup g 2 =1 K 1 -e -h (g), g = sup g 2 =1 lim p→+∞ K 1 -e -hp (g), g ≤ lim p→+∞ sup g 2 =1 K 1 -e -hp (g), g = lim p→+∞ sup g 2 =1 K 1 -e -hp g , 1 -e -hp g Because √ 1 -e -hp g has a compact support, K 1 -e -hp g , 1 -e -hp g = K |B(0,p) 1 -e -hp g , 1 -e -hp g ,
where K |B(0,p) is the restriction of K on L 2 (B(0, p)). Since K |B(0,p) is trace-class, we have

K |B(0,p) 1 -e -hp g , 1 -e -hp g ≤ λ max p 1 -e -hp g 2 ≤ λ max p ≤ λ max ,
where λ max p (resp. λ max ) is the greatest eigenvalue of K |B(0,p) (resp K). Then sup

g 2 =1 K 1 -e -hp g , 1 -e -hp g ≤ λ max and lim p→+∞ sup g 2 K 1 -e -hp g , 1 -e -hp g ≤ λ max . Thus, K 1 -e -h ≤ λ max < 1
which proves the result.

So we have an upper bound of 1 n

Tr K 1 -e -hp n independent of p which is summable so we can exchange the limit and the sum.

Secondly, we prove that lim

p→+∞ Tr K 1 -e -hp n = Tr K 1 -e -h n .
For all n ≥ 0:

Tr K 1 -e -h n -Tr K 1 -e -hp n = Tr K 1 -e -h n -K 1 -e -hp n = Tr K 1 -e -h -K 1 -e -hp n-1 k=0 K 1 -e -h k K 1 -e -hp n-1-k ≤ n-1 k=0 Tr K 1 -e -h -K 1 -e -hp K 1 -e -h k K 1 -e -hp n-1-k .

Moreover we have the following inequalities

Tr K 1 -e -h -K 1 -e -hp K 1 -e -h k K 1 -e -hp n-1-k ≤ Tr K 1 -e -h -K 1 -e -hp K 1 -e -h k K 1 -e -hp n-1-k ≤ Tr K 1 -e -h -K 1 -e -hp K 1 -e -h k K 1 -e -hp n-1-k ≤ Tr K 1 -e -h -K 1 -e -hp K 1 -e -h n-1 ,
the last inequality taking place according to (4.1). We finally have the following inequality

Tr K 1 -e -h n -Tr K 1 -e -hp n ≤ Tr K 1 -e -h -K 1 -e -hp n K 1 -e -h n-1 -→ p→+∞ 0. lim p→+∞ exp - ∞ n=1 1 n Tr K 1 -e -hp n = exp -lim p→+∞ ∞ n=1 1 n Tr K 1 -e -hp n = exp - ∞ n=1 1 n lim p→+∞ Tr K 1 -e -hp n = exp - ∞ n=1 1 n Tr K 1 -e -h n .
In order to prove the convergence in distribution of n(ρ) -1 M ρ (µ), for µ ∈ M + α,β , we study the convergence of its Laplace transform : for θ ≥ 0,

E exp -θn(ρ) -1 M ρ (µ) = exp θE[n(ρ) -1 M ρ (µ)] E exp -θn(ρ) -1 M ρ (µ) = exp θE[n(ρ) -1 M ρ (µ)] E exp - R d ×R + ×R + θn(ρ) -1 mµ(B(x, r))Φ ρ (dx, dr, dm) .
To compute this last term, we use Proposition 3.1. The hypothesis (11) in this proposition is satisfied because in our context of weighted balls model, we have h(x, r, m) = mµ(B(x, r)) and therefore:

R d ×R + ×R + (1-e -mµ(B(x,r)) )K ρ (x, x)f (r/ρ)g(m)dx dr ρ dm ≤ λ(ρ) sup x∈R d K(x, x) R d ×R + ×R + mµ(B(x, r))f (r/ρ)g(m)dx dr ρ dm ≤ λ(ρ)ρ d v d µ(R d ) sup x∈R d K(x, x) R + mg(m)dm R + r d f (r)dr < +∞.
Moreover, we have to check that

K ρ 1 -e -h ∈ L 2 (R d × R + × R + )
, if we denote by h the function given by h(x, r, m) mµ(B(x, r)) defined on

R d × R + × R + . (R d ×R + ×R + ) 2 K ρ 1 -e -h 2 ((x, r, m), (y, s, m ′ ))dxdrdmdydsdm ′ = (R d ×R + ×R + ) 2 1 -e -mµ(B(x,r)) g(m) f (r/ρ) ρ K 2 ρ (x, y) × f (s/ρ) ρ g(m ′ ) 1 -e -m ′ µ(B(y,s)) dxdrdmdydsdm ′ ≤ (a) λ(ρ) 2 sup x∈R d K(x, x) 2 R d ×R + ×R + 1 -e -mµ(B(x,r)) g(m) f (r/ρ) ρ dxdrdm 2 ≤ λ(ρ) 2 sup x∈R d K(x, x) 2 ρ d v d µ(R d ) R + mg(m)dm R + r d f (r)dr 2 < +∞,
where inequality (a) stands because

K 2 ρ (x, y) ≤ K ρ (x, x)K ρ (y, y) ≤ λ(ρ) 2 sup x∈R d K(x, x) 2
thanks to [START_REF] Kaj | Scaling limits for random fields with long-range dependence[END_REF].

We can now apply Proposition 3.1 to obtain:

E exp - R d ×R + ×R + θn(ρ) -1 mµ(B(x, r))Φ ρ (dx, dr, dm) = exp - ∞ n=1 1 n Tr K 1 -e -θn(ρ) -1 h n . (17)
The Laplace transform of n(ρ) -1 M ρ (µ) thus rewrites E e -θn(ρ) -1 Mρ(µ) = exp

R d ×R + ×R + ψ θn(ρ) -1 mµ(B(x, r)) λ(ρ)K(0) f (r/ρ) ρ g(m)dxdrdm × exp - n≥2 1 n Tr K ρ 1 -e -θn(ρ) -1 h n (18) with ψ(u) = e -u -1 + u.
The convergence of (18) derives from the following lemmas. The complete proofs of these lemmas are given in [3].

Lemma 3.2 For all n ≥ 2, we have

Tr K ρ 1 -e -θn(ρ) -1 h n ≤ Tr K ρ 1 -e -θn(ρ) -1 h 2 n/2 .

Proof:

The key point is that K ρ 1 -e -θn(ρ) -1 h is trace-class, and as a consequence a Hilbert-Schmidt operator thanks to Lemma 4.2 so the proof of Lemma 2.9 in [3] applies in the same way.

Assume Conditions

(2), ( 6) and (8), and consider µ ∈ M + α,β . Then there is a constant M ∈]0, +∞[ such that, Note that, in this non-stationary setting, the Poissonian limits for n = 1 come from Theorem 1, Theorem 2 and Theorem 3 in [START_REF] Gobard | Random balls model with dependence[END_REF] taking f (x, r) = K(x, x)f (r) in our situation.

Tr K ρ 1 -e -θn(ρ) -1 h 2 ≤ Mθ 2 λ(ρ)ρ q n(ρ)
As in [3], it is now enough to show now that lim ρ→0

Tr K ρ 1 -e -θn(ρ) -1 h 2 = 0 in the three regimes.

As a consequence of Lemma 3.3, it remains to show that λ(ρ)ρ q n(ρ) 2 -→ ρ→0 0.

(i) Large-balls scaling. Since lim ρ→0 λ(ρ)ρ β = +∞, for ρ small enough we have λ(ρ)ρ β ≥ 1 and so λ(ρ)ρ β 1/α ≥ λ(ρ)ρ β 1/2 with α ∈]1, 2]. Thus since q > β we have : 0 ≤ λ(ρ)ρ q n(ρ) 2 ≤ λ(ρ)ρ q λ(ρ)ρ β = ρ q-β -→ ρ→0 0.

(ii) Intermediate scaling. In this case, n(ρ) = 1 and since q > β we have : 0 ≤ λ(ρ)ρ q n(ρ) 2 = λ(ρ)ρ q = λ(ρ)ρ β ρ q-β -→ ρ→0 0.

(iii) Small-balls scaling. Since we consider µ ∈ L 1 (R d ) ∩ L 2 (R d ), Proposition 2.5-(ii) in [3] (which remains correct for µ ∈ L 1 (R d ) ∩ L 2 (R d )) ensures that we can take q = 2d and then with n(ρ) = (λ(ρ)ρ β ) 1/γ since β < αd ≤ 2d we have : 0 ≤ λ(ρ)ρ 2d n(ρ) 2 = λ(ρ) (β-2d)/β -→ ρ→0 0.

4 Appendix: Lemmas for the Laplace of DPP

In this appendix, we state and prove three different lemmas used in section 3 to prove Proposition 3.1.

Lemma 4.1 If f, g are two real functions on E such that 0 ≤ f ≤ g, then we have

K [f ] ≤ K [g]
where K [f ] is the operator with kernel K [f ] (x, y) = f (x)K(x, y) f (y) for x, y ∈ E.

Proof: Recall that:

K [f ] = sup h∈L 2 (E) K [f ] (h), h h 2 2 .
Let h ∈ L 2 (E). where l = f /g ≤ 1. So we have:

K [f ] (h), h = E K [f ] (h)(x)h(x)dx = E E f (
K [f ] (h), h ≤ K [g] lh 2 2 ≤ K [g] h 2 2
and the result follows.

We recall the following result from Proposition 280 in [7]: Lemma 4.2 A trace-class operator K on a Hilbert space H is a compact operator.

Proposition 1 . 1

 11 Assume (5) is in force. For all µ ∈ Z(R d ), E [M(|µ|)] < +∞. As a consequence, M(µ) in (4) is almost surely well defined for all µ ∈ Z(R d ).

  x)K(x, y) f (y)h(y)h(x)dydx = )h(x)dydx = K [g] (lh) , lh

  2 . We give a short proof of Theorem 2.2. In this non-stationary case, the proof follows the same general strategy as in[3]. Roughly speaking, the limits are driven by the term n = 1 in (17) while the other terms (n ≥ 2) are still negligible.

	Proof: The computations are analogous to that in [3] for the model without weight. It is
	important to observe that the key point, namely inequality (25) in [3], remains true because
	µ ∈ M + α,β so thanks Proposition 2.2 (iii) in [4], µ ∈ M + 2,β = M + β using the notations of [3]. 2
	To be complete, the constant M is equal to C K C µ C f	mg(m)dm	with the notations
	R +
	of Lemma 2.10 in [3].
	Proof:


Cristian Felipe Ocampo-Blandon

Yann Gousseau
email: yann.gousseau@telecom-paristech.fr

Non Local Exposure Fusion

Keywords: exposure fusion, image fusion, patches, non-local methods, high dynamic range

come

Introduction and Previous Works

Natural scenes often exhibit luminance ranges that are beyond the capacity of most imaging sensors. In such cases, a single standard photograph cannot capture details in dark regions without producing saturated values in the brightest regions. A classical way to create a high dynamic range (HDR) image faithfully representing the luminance of a scene is to combine photographs acquired with different exposure times [START_REF] Debevec | Recovering high dynamic range radiance maps from photographs[END_REF]. Once the response function of the camera is known (or when working with linear RAW images) and when the scene is static, this combination task is essentially a statistical problem [START_REF] Aguerrebere | Best algorithms for HDR image generation. a study of performance bounds[END_REF]. In order to visualize such an image on a classical display, contrast is then compressed by using a so-called tone mapping operator [START_REF] Fattal | Gradient domain high dynamic range compression[END_REF][START_REF] Durand | Fast bilateral filtering for the display of highdynamic-range images[END_REF]. An alternative to this two-step procedure (HDR image creation followed by tone mapping) was proposed in [START_REF] Mertens | Exposure fusion[END_REF] and is called exposure fusion. The main idea is to bypass the HDR creation step and to directly create a low dynamic range image (typically made of 8 bits per color channel) by fusing the input images corresponding to different acquisition times. Specific weights are used to ensure that the final result is contrasted enough, has vivid colors and avoid over and under-exposure. Using classical ideas from computer graphics [START_REF] Ogden | Pyramid-based computer graphics[END_REF], the images are fused in a multi-scale framework, enabling one to blend images seamlessly.

This procedure is very efficient and has yielded numerous softwares and plugins. It has also triggered an abundant research literature proposing variants on the original method: involving image decomposition into several components at the patch level [START_REF] Ma | Multi-exposure image fusion: A patch-wise approach[END_REF], formulating the problem as a MAP estimation [START_REF] Song | Probabilistic exposure fusion[END_REF] or involving the bilateral filter [START_REF] Raman | Bilateral filter based compositing for variable exposure photography[END_REF]. Nevertheless, such approaches have a common drawback: they fail when there is either camera shake or moving objects in the scene. Indeed, they all perform the fusion at pixel level and assume that images are registered and that objects are still. Of course, this is a serious limitation in practice, and several works have tackled this issue. Anti-ghosting algorithms [START_REF] An | A multi-exposure image fusion algorithm without ghost effect[END_REF][START_REF] Pece | Bitmap movement detection: HDR for dynamic scenes[END_REF][START_REF] Zheng | Hybrid patching for a sequence of differently exposed images with moving objects[END_REF] propose to perform image alignment and possibly to explicitly detect moving objects to prevent using them in the fusion. Recently, the methods from [START_REF] Hu | HDR deghosting: How to deal with saturation?[END_REF][START_REF] Qin | Robust match fusion using optimization[END_REF] proposed to solve the problem by an iterative optimization procedure relying on correspondences between patches. In particular, the authors of [START_REF] Hu | HDR deghosting: How to deal with saturation?[END_REF] propose to create, from a reference image and a set of differently exposed images, a set of images that are all aligned and radiometrically coherent with the reference, and whose content is automatically modified (in the case of moving objects) to match the reference. This is achieved thanks to contrast prescription and patch-based content comparison between the images. As illustrated in [START_REF] Hu | HDR deghosting: How to deal with saturation?[END_REF], classical exposure fusion can then be applied to the aligned set to yield a satisfactory final image, even in the case of camera shake and object motion.

In this paper, we propose an exposure fusion method that also deals both with camera shake and object motion. The basic idea is very simple and illustrated by Figure 1: instead of fusing values taken by pixels at the same spatial position, we fuse values of pixels having similar neighborhood (thereafter called a patch) in the different images. The method is therefore in the spirit of non-local restoration methods such as the non-local means [START_REF] Buades | A non-local algorithm for image denoising[END_REF], and also share similarities with the non-local method introduced in [START_REF] Aguerrebere | Simultaneous HDR image reconstruction and denoising for dynamic scenes[END_REF] for the creation of HDR images. Such a simple approach has the ability to deal with camera motion (although a previous global image registration may be useful to accelerate the process or to enhance the similarity between patches) and with moving objects. By construction, the method also boils down to the original exposure fusion algorithm in the absence of camera and object motions. The rest of the paper is organized as follows. In Section 2, the original exposure fusion algorithm is first briefly recalled and a detailed presentation of the proposed method is provided. In Section 3, we show experimental results and comparisons with the recent state-of-the-art algorithm from [START_REF] Hu | HDR deghosting: How to deal with saturation?[END_REF]. 2 Non-local Exposure Fusion

Classical Exposure Fusion

The method introduced in [START_REF] Mertens | Exposure fusion[END_REF] proposes to fuse a series of images I 1 , . . . , I N acquired with different exposure settings (the knowledge of these is not needed).

For each pixel x and index i a weight W i (x) is defined by taking into account the quality of contrast, color saturation and well-exposedness at this pixel. The idea is then to fuse the values I i (x) according to the weights. The most straightforward approach would be to define the resulting image R as

R(x) = N i=1 W i (x)I i (x) , (1)
but this yields incoherences in flat regions and visible seams at slow transitions.

In order to achieve seamless fusion, the blending is performed in a multi-scale framework. For each level l of a Laplacian pyramid [START_REF] Ogden | Pyramid-based computer graphics[END_REF], the Laplacian pyramid of the resulting image R is computed as

L l (R) = N i=1 G l (W i )L l (I i ) ,
where L l (I i ) is the Laplacian pyramid at level l of image i and G l (W i ) is the Gaussian pyramid at level l of the weight map W i . The final image R is then reconstructed from its Laplacian pyramid. As explained in the introduction, this (otherwise very efficient) method fails when there is camera shake or object motions.

Non-local Fusion

In order to deal with dynamic scenes we propose, for each pixel x, to replace the fusion of values at this pixel with the fusion of values at well chosen pixels f 1 (x), . . . , f N (x) in the series of images. In the static case these pixels would simply be all equal to x, and in the general case, they will be defined in order to compensate motions, as we will now see. In a way similar to the HDR creation method from [START_REF] Aguerrebere | Simultaneous HDR image reconstruction and denoising for dynamic scenes[END_REF], we define f i (x) to be the pixel having its neighborhood in I i the closest to the one of x in I ref (where ref ∈ (1, . . . , N )) a reference image to be defined later. More precisely, we define P i (x), a patch at pixel x, to be the collection of values in I i in a (2s + 1) × (2s + 1) square window centered at x. We assume that we have a distance d between patches (to be defined in the next section) and then define the pixel most similar to x in I i to be

f i (x) = arg min y∈x+N (d(P ref (x), P i (y))) , (2)
where N is a search window (bigger than the patch). In the (unlikely) case where the minimum is achieved for several values, one is chosen at random. If the minimum is achieved for only one value, we have in particular that f ref (x) = x.

Let us now assume for a moment that we perform the fusion at a single scale (that is, as in Formula (1)). Then the resulting image could be defined as

R(x) = N i=1 W i (f i (x))I i (f i (x)) •
Of course, as in the classical exposure fusion, defining the resulting image R this way would yield incoherences and seams. Therefore, the non-local fusion is defined, for each level l of a multi-scale pyramid, as:

L l (R) = N i=1 G l ( Wi )L l ( Ĩi ) , (3)
where Ĩi is defined as Ĩi (x) = I i (f i (x)). The weights Wi are obtained by multiplying the original weights w 1 i from [START_REF] Mertens | Exposure fusion[END_REF] (accounting for well-exposedness, color saturation and contrast) at position f i (x) by a weight asserting the similarity between x and f i (x). Similarly to what is done in the NL-means algorithm [START_REF] Buades | A non-local algorithm for image denoising[END_REF] we define this similarity weight to be

w 2 i (x) = exp - d(P ref (x), P i (f i (x))) 2 σ 2 , ( 4
)
where σ is a parameter. The resulting weights are Wi (x) = w 1 i (f i (x)) × w 2 i (x). In short, for a given pixel, the fusion is performed by using the most similar pixels in each image of the sequence and by modulating the original weights of the fusion by a similarity score.

Computing Pixel Similarity

Here, we detail how to compute the distance between patches. For this, images acquired with different exposure settings (e.g. different exposure times) should be comparable. In this paper, we assume that we have access to the linear RAW images of the captured scene. This is not a real limitation in the common case of embedded processing. When these RAW images are not available, other scenarios are possible (such as matching all image histograms to the one of the reference image, in a way similar to [START_REF] Hu | HDR deghosting: How to deal with saturation?[END_REF]), but these will not be investigated in this paper.

In addition to the processed images I 1 , . . . , I N (typically JPEG images after demosaicking, white balance and gamma transform), we consider the corresponding RAW (linear) images R 1 , . . . , R N . We assume that these have been acquired with exposure times t 1 , . . . , t N . The luminance images can then be computed as

L i = gt -1 i (R i -O),
where O is the offset of the camera and g the gain (see e.g. [START_REF] Aguerrebere | Study of the camera acquisition process and statistical modeling of the sensor raw data[END_REF]). Without loss of generality, we will assume that g = 1 and that the offset is known. The comparison of pixels is then carried out in the luminance domain, so that the distance d between two patches P i (x) and P j (y) is defined as the Sum of Squared Differences (SSD) between the values of the patches in the respective luminance images L i and L j , that is, If the camera offset is unknown, or in order to achieve results that are more robust to spatial variability of the gain (photo response non-uniformity, PRNU) or imprecision in the exposure time, the L 2 distance can be replaced by a similarity measure that is invariant to affine transforms of the luminance. Although this aspect will not be developed here, we have observed that the affine invariant distance proposed in [START_REF] Delon | Stabilization of flicker like effects in image sequences through local contrast correction[END_REF] is very robust to errors in the luminance conversion.

d(P i (x), P j (y)) = s k=-s (L i (x + k) -L j (y + k)) 2 • ( 5
) (a) (b)

Algorithm Details

First, the images are aligned through a global SIFT-based homography estimation and bicubic interpolation. Although not necessary, we observed that this step yields more reliable patch comparisons. A reference image is picked as in [START_REF] Hu | HDR deghosting: How to deal with saturation?[END_REF] by choosing the image having the least number of saturated pixels. Then, saturated regions in the reference are removed (see below). Second, and for each pixel in the reference, its best match in image i is chosen using Formula (2) and the similarity measure defined by [START_REF] An | A multi-exposure image fusion algorithm without ghost effect[END_REF]. In order to accelerate the algorithm, the best match is found using the PatchMatch algorithm [START_REF] Barnes | Patchmatch: A randomized correspondence algorithm for structural image editing[END_REF], yielding, for each patch of the reference image, an approximate nearest neighbor (ANN) in each image L i . The PatchMatch algorithm is run with its default parameters, which in particular implies that the search window N is the entire image.

The corresponding value of the final image is obtained by fusing the values of the best matches using Formula (3). In this formula, the weights Wi are obtained by multiplying the original weight of [START_REF] Mertens | Exposure fusion[END_REF] by the similarity weight defined by Formula (4). The original weights from [START_REF] Mertens | Exposure fusion[END_REF] depend on three exponents that we all choose to be equal to 1. The parameter σ of the similarity measure (4) is chosen as the quantile of level 5% of the distribution of the distances between patches (separately for each couple L ref , L i ). We empirically found this choice to provide a reasonable compromise in order to simultaneously obtain a result similar to the original algorithm from [START_REF] Mertens | Exposure fusion[END_REF] and efficiently discard moving objects.

A last point that we did not explain is what to do with saturated pixels. We could have chosen a strategy similar to the one of [START_REF] Hu | HDR deghosting: How to deal with saturation?[END_REF]: use only the non-saturated pixels when comparing patches. However, this fails as soon as the saturated region is wider than the patch, which is very common in practice. In order to get more robust results, we use the same strategy as in [START_REF] Aguerrebere | Simultaneous HDR image reconstruction and denoising for dynamic scenes[END_REF]. For each connected (e) Our Method -NLEF, (184 sec). (f) Fusion with [START_REF] Hu | HDR deghosting: How to deal with saturation?[END_REF], (394 sec). region Ω of saturated pixels, we detect the region from the other images that is the most similar to Ω on its immediate neighborhood (a ring outside Ω). This region is then blended into the reference using Poisson editing [START_REF] Pérez | Poisson image editing[END_REF].

Experiments

In all experiments (see Figure 3), we use a search window equal to the full image. The patch size is set to s = 1 (that is 3 × 3 patches).

First, as a sanity check, we compare the result of our method, NLEF (nonlocal exposure fusion), with that of the original exposure fusion [START_REF] Mertens | Exposure fusion[END_REF] on a static scene with a moving object, from the database [START_REF]international university of sarajevo[END_REF]. As we can see in Figure 2, the results are very similar on static parts (everything except one object on the table), but no ghost is produced with our method for the moving object.

Then, we provide three comparisons between our method and the recent method from [START_REF] Hu | HDR deghosting: How to deal with saturation?[END_REF], using the code kindly provided by the authors. Since exposure fusion does not aim at producing a real HDR image (a faithful representation of the luminance values), we only provide visual comparisons. A full user study would be the best way to evaluate the quality of the results [START_REF] Ma | Perceptual quality assessment for multiexposure image fusion[END_REF] but is beyond the scope of this paper. Both methods rely on PatchMatch to get patch correspondences, and we use the same parameters for this algorithm in all experiments and for the two approaches (full search over the complete image and patch size of s = 1). The examples corresponding to Figures 3(a) and 3(c), are acquired with a Canon 400D and Canon 7D camera, respectively, for which the offsets are 256 and 2046. Images of these experiments, as well as other comparisons, are included on the supplementary material.

In Figures 3(a) and 3(e), one can see that our method produces more accurate colors. However, the color errors from [START_REF] Hu | HDR deghosting: How to deal with saturation?[END_REF] could probably be solved by taking into account the luminance computed from linear images. Second, one observes that both methods produce accurate details and no ghosts, with occasional artifacts at different positions. However, our approach is both faster and much simpler, and in particular involves no complex optimization procedure.

Conclusions

In this work, we have introduced a non-local exposure fusion yielding good results in the presence of camera shake and objects' motions. There are several ways this work could be continued. First, the saturated parts are handled using a Poisson editing procedure, but this step could probably be integrated in the PatchMatch algorithm, at the propagation step. Second, the method may produce noisy results, in particular when the reference image is noisy. We plan to tackle this issue by using more patches for the fusion, in a way similar to the NL-means algorithm. Last, similar non-local approaches could be studied in different fusion schemes, for instance to perform focus stacking.

Fig. 1 .

 1 Fig. 1. Instead of fusing values at the same position (a), we use a non-local fusion (b) with values from different positions, to account for object motions and camera shake.

Fig. 2 .

 2 Fig. 2. Image fusion: (a) non-local exposure fusion (NLEF), (b) exposure fusion [15].

  (a) Our Method -NLEF, (69 sec).(b) Fusion with[START_REF] Hu | HDR deghosting: How to deal with saturation?[END_REF], (143 sec).(c) Our Method -NLEF, (84 sec). (d) Fusion with[START_REF] Hu | HDR deghosting: How to deal with saturation?[END_REF], (116 sec).

Fig. 3 .

 3 Fig. 3. Sets of bracketed exposure images and their fusion. (Third image set is courtesy of EMPA Media Technology.)

Acknowledgments. The authors thank C. Aguerrebere for her help and advice. This work is supported by Colciencias under the grant "Programa doctoral de Becas" call 529 and the French Research Agency (ANR) under project MIRIAM.


H Lavaste

C Pillet

Sauret

S Hybois

A Lombart

P Puchaud

J Bascou

F Lavaste

H Pillet

C Sauret

Effects of ellipsoid parameters on scapula motion during manual wheelchair propulsion based on multibody kinematics optimization. A preliminary study

Keywords: musculoskeletal, shoulder, scapulothoracic, wheelchair, kinematics

Effects of ellipsoid parameters on scapula motion during manual wheelchair propulsion based on multibody kinematics optimization. A preliminary study

Introduction

Recording the motion of the scapula during manual wheelchair (MWC) propulsion is of particular interest but remains a challenge due to the soft tissue artefact. For that purpose, several authors recommend the use of multibody kinematics optimization (MKO) [START_REF] Naaim | Effect of Various Upper Limb Multibody Models on Soft Tissue Artefact Correction: A Case Study[END_REF] but the way to adapt the generic model to the subject (scaling) is not as trivial as it can appear [START_REF] Duprey | Kinematic Models of the Upper Limb Joints for Multibody Kinematics Optimisation: An Overview[END_REF]. In the description of the upper limbs, the main challenge is to accurately model the shoulder complex, including sternoclavicular, scapulothoracic, acromioclavicular and glenohumeral joints. In the field of MWC propulsion research, several upper extremity models have been used. The simplest had modelled the shoulder complex as a ball joint between the thorax and the humerus located on the glenoid centre. A more complex model implemented regression equations allowing sternoclavicular and scapulothoracic coupled coordinates [START_REF] Holzbaur | A Model of the Upper Extremity for Simulating Musculoskeletal Surgery and Analyzing Neuromuscular Control[END_REF]. Finally, others authors had proposed to guide the motion of the scapula on an ellipsoid to represent the scapulothoracic joint [START_REF] Van Der Helm | A Finite Element Musculoskeletal Model of the Shoulder Mechanism[END_REF][START_REF] Seth | A Biomechanical Model of the Scapulothoracic Joint to Accurately Capture Scapular Kinematics during Shoulder Movements[END_REF]. However, the personalization of the ellipsoid parameters is not obvious.

The aim of this work was to evaluate the influence of the ellipsoid parameters (centre location and radii) on kinematics reconstructed using MKO during MWC propulsion.

Methods

Experimental data

One able-bodied subject (22 y.o, 1.75m, 62 kg) participated in this preliminary study. During the experiment, the subject was asked to perform several propulsion cycles on a MWC ergometer. 36 reflective markers were placed on the upper limbs, and their 3D locations were recorded at 100 Hz by an 8-cameras optoelectronic motion capture system (Vicon system, hardware: T40 cameras; Nexus software; Oxford Metrics, UK).

Musculoskeletal model

A three-dimensional linked-segment model of the upper limbs was used to compute biomechanical parameters. The model included hand, forearm bones, humerus, clavicle, scapula and thorax. The latter were derived from the generic scapulothoracic joint model [START_REF] Seth | A Biomechanical Model of the Scapulothoracic Joint to Accurately Capture Scapular Kinematics during Shoulder Movements[END_REF], which was symmetrised using a custom routine. This model was scaled to the subject using the algorithm implemented in OpenSim software [START_REF] Delp | OpenSim: Open-Source Software to Create and Analyze Dynamic Simulations of Movement[END_REF], and the location of the acromio-clavicular contact point was scaled with a custom program.

Data processing

Joint kinematics was obtained for a steady-state propulsion cycle, using an inverse kinematics algorithm implemented in OpenSim software. The impact of the estimation of scapulothoracic joint parameters was evaluated with a full factorial design with 6 parameters on two levels. Parameters were the 3 radii of the ellipsoid and the 3 coordinates of its centre in the thorax frame; and were alternately defined at 80% and 120% of the initial values obtained from the default scaling step.

For each of the 64 resulting models, the overall root mean squared error (RMSE) between reconstructed and experimental markers was computed, and compared to the value obtained with the initial default scaled model.

Results and discussion

The results showed that among the 64 models, only eight induced a lower error between the reconstructed and the experimental position of the scapula markers with respect to the default scaled model. As depicted in Figure 1, most of the modified models induced a lower RMSE during a propulsion cycle for the arm (number of such case: 34), forearm (33) and hand (51), reaching up to a 5% improved reconstruction for the arm. A lower number of markers

Conclusions

A full factorial design showed that scapular ellipsoid parameters can result in significant improvement in the kinematics reconstruction during MWC propulsion. In this perspective, subject-specific imaging methods coupled to optimization algorithms could be of great interest. This study could also be extended with the clavicle length which was already showed as critical [START_REF] Duprey | Kinematic Models of the Upper Limb Joints for Multibody Kinematics Optimisation: An Overview[END_REF].

improved the kinematics reconstruction of the thorax ( 12) and clavicle ( 25), but the clavicle markers reconstruction reached 9% of improvement depending on the ellipsoid parameters chosen. Among the combinations that enhanced the scapula markers, several ellipsoid parameters led to a better RMSE for multiple body segments, with two distinct situations: either improving thorax/scapula only, or improving all the segments except the thorax.

The impact of ellipsoid parameters on generalized coordinates is shown on Figure 2, for two selected degrees of freedom of the scapula, namely abduction and elevation. The ellipsoid parameters from the OpenSim default scaling (plain blue line) were obviously within the standard deviation of the 64 modified parameters (light blue space), whereas the parameters that optimized the global marker RMSE (red dashed line) displayed a different behaviour along the propulsion cycle.

Since many models resulted in higher RMSE than the default model (Figure 1), it could be assumed that the default scaling method implemented in OpenSim provided a reasonable kinematics reconstruction. However, by slightly changing ellipsoid parameters, it was shown that RMSE could be reduced. With the hypothesis that a lower marker error denotes a better kinematic chain, this result implied that ellipsoid parameters are key elements that should be subject-specifically defined.

The differences in kinematics can be explained from geometric interpretation of ellipsoid parameters: the ellipsoid that minimized the global marker RMSE placed the scapulothoracic joint centre more downwards and backwards in the thorax, which explained why the resulting scapula elevation is higher (Figure 2). Moreover, for some degrees of freedom, like scapula abduction, the standard deviation resulting from the different models was reduced around the middle of the propulsion cycle, whereas it increased at the beginning or the end of the propulsion cycle. This could be explained by the scapula motion that exhibit a lower amplitude in the middle of the propulsion cycle, while it reaches greater angles when the arm moves backwards to allow the rim grasp.

Figure 1 .

 1 Figure 1. distribution of the ratio between the markers rmse obtained for the 64 models and with the default scaled model. results were presented for each segment. Value lower than 1 means improvement through the modified model

Figure 2 .

 2 Figure 2. impact of ellipsoid parameters on scapula angles during mWC propulsion


J Bascou

C Fallot

H Pillet

S Hybois

F Lavaste

C Sauret

Assessment of power losses due to ground contact forces during usual manual wheelchair movements

Keywords: manual wheelchair, movement resistances, rolling, turning

Introduction

Manual wheelchair (MWC) allows disabled people to recover autonomy, but may overload their upper limbs, causing pain or musculoskeletal disorders [START_REF] Mercer | Shoulder joint kinetics and pathology in manual wheelchair users[END_REF]. The stresses sustained by the MWC user are due to the muscular forces necessary to execute useful MWC movements (accelerating in straight line or in rotation, keeping up the velocity, climbing up a ramp …). Unfortunately, a non-negligible part of these efforts is wasted in useless movement resistances [START_REF] Lin | Evaluation of wheelchair resistive forces during straight and turning trajectories across different wheelchair configurations using freewheeling coast-down test[END_REF]. Thus, lowering these resistances, mainly due to ground/wheel contact, would favour the decrease of the constrains sustained by the user upper limbs. Many studies focused on assessing these resistances for straightforward propulsion of the MWC, but very few addressed the quantification of the energy lost during turning manoeuvers [START_REF] Lin | Evaluation of wheelchair resistive forces during straight and turning trajectories across different wheelchair configurations using freewheeling coast-down test[END_REF]. Recently, our team quantified the values of wheel swivelling resistance torque, which depends on the wheel type, ground nature, load applied on the wheel and the curvature radius of the wheel trajectory (Fallot et al. under review). Based on these values, the aim of the present study was to quantify the energy lost during typical daily living activities with a MWC.

Methods

Subjects and test protocol

Eight healthy subjects propelled an instrumented MWC (FRET 2), equipped with 6-components dynamometers on both handrims and on the chassis, rear wheel angular potentiometers and an inertial measurement unit (Sauret et al. 2012) during 5 min blocks of level-ground displacement (covered with low-pile carpet), performed at self-selected speed. Two measurement sessions of three blocks were performed by all subjects. Each block consisted in a guided succession of start-ups, propulsions in straight line, slaloms, stops, and half turns.

Data gathering and analysis

Using the instrumented MWC, the user efforts on both handrim, the seat, the backrest and the footrest, as well as the MWC kinematics could be measured. Vertical ground reaction forces on front and rear wheels were assessed as described in [START_REF] Sauret | Proposal of an index for evaluating pitch instability during actual locomotion with a manual wheelchair[END_REF]. Each wheel kinematics were also assessed. The mass and center of mass fore-aft position of the instrumented MWC was previoulsy determined using force plateforms.

By considering each wheel type and radius, its rolling resistance factor could be assessed (Sauret et al. 2012). By considering each wheel type and MWC instantaneous center of rotation (ICR), its swivelling resistance factor could be determined. Multiplying rolling and swivelling resistance factors by the load applied on each wheel gave the rolling and swivelling resistances torques. By multiplying these torques by each wheel angular velocity, the power losses could be computed. Time integration of the power losses provided the energy spent during each activity.

MWC movements were classified in straightforward propulsion (linear velocity superior to 0.02 m/s and rotation velocity inferior to 0.2 rad/s), and turning motion (rotation velocity superior to 0.2 rad/s).

The mean power lost in rolling and swivelling resistances, for each wheel and in each movement classification, were then compared.

Results and discussion

The instrumented MWC mass was 38.2 kg and its centre of mass fore-aft position was 10.9 cm from the rear axle. 100 N load applied on front or rear wheels induced rolling resistance torques of 3.5 and 1.5 Nm respectively. Furthermore, considering in addition a 0.4 m gyration radius for the wheel trajectory, front and rear wheel swivelling resistance torques were 0.91 and 0.084 Nm respectively. similarly distributed between front and rear wheels, high discrepancies were found between wheels during MWC rotation. This prevalent effect of rear wheel in rotation movement could be interesting when setting up a MWC, as front wheels could be chosen only considering rolling resistance performance. As rear wheels induce energy losses in both manoeuver types, their rolling and swivelling performances should be considered.

Healthy subject movements certainly differed from those of MWC experts, but for comfortable speed, the effect of this difference on front and rear ground reaction force and speed may not be enough to alter significantly the conclusions of this study.

Computing the mean energy spent in 1 min trial in rolling and swivelling resistances, which could be roughly compared to 1 min of daily life activity, showed that the user had to produce 1112 J to overcome the main MWC movement resistances.

Conclusions

The knowledge individual wheel contact forces between MWC and ground was used to assess power losses in any situation. It also opens the path to realistic simulation of MWC propulsion, not limited to straightforward movements. Power losses in MWC turning motion proved to be non-negligible: as they may often occur during daily living (propulsion in an apartment, sports …), they should be addressed when trying to adapt the MWC settings to the user environment.

Averaged ground reaction forces on front / rear wheels (all subjects and trials mixed up) were respectively 281 N/817 N for straightforward movement and 291 N/797 N for rotation. These values were similar for both movements and could be explained by the choice of comfortable self-selected speed which imposed few trunk movements resulting in few weight transfers on front wheels. The larger weight applied on front wheels during rotation movement could be explained by the bending of the trunk during this movement. The total time spent in linear movements (275 s, summing all trial and subject durations) and to the one spent in rotation manoeuvers (218 s) were similar.

The mean of linear velocities performed were 0.70 m/s for straightforward movement and 0.78 m/s for rotation movements, which means that most of the rotation movements were not 'on-the-spot' rotations. Mean angular velocity was 0.18 rad/s (10°/s) for rotation movements. These values were inferior to the velocity usually observed during straightforward propulsion as start-up phases were considered in the movements, where initial velocity was near to zero.

Mean power losses during straightforward movement were 7.6 W by front wheels and 9.4 W by rear wheels. Indeed, the rear wheel rolling factor was lower than the front wheel one but higher load was applied on rear wheels, increasing the power they dispelled (Sauret et al. 2012).

Concerning MWC rotation manoeuvers, the mean of power losses due to wheel swivelling resistance were 0.1 W for front wheels and 4.3 W for rear wheels. These results were not surprising considering that the front wheel swivelling resistance parameter was far lower than the rear wheel one and that the load applied on front wheels was also lower than the one applied on rear wheels. During turning manoeuvers, power losses due to rolling resistance in front and rear wheels were respectively 5.8 and 6 W, which is comparable to power losses in wheel swivelling.

Swivelling resistance proved to be non-negligible with respect to rolling resistance and should be considered when quantifying the energetic performance of a MWC. Whereas rolling resistance power losses proved to be

Figure 1 .

 1 Figure 1. mean power of rolling and swivelling resistances, for each movement and wheel types.

Acknowledgements

The authors would like to thank the 'Fondation du sport français -Henri Sérandour' for their financial support.