Oudom Somphone
email: oudom.somphone@philips.com

Benoit Mory
email: benoit.mory@philips.com

Sherif Makram-Ebeid
email: sherif.makram-ebeid@philips.com

Laurent Cohen
email: cohen@ceremade.dauphine.fr

Prior-based Piecewise-smooth Segmentation by Template Competitive Deformation using Partitions of Unity

We propose a new algorithm for two-phase, piecewise-smooth segmentation with shape prior. The image is segmented by a binary template that is deformed by a regular geometric transformation. The choice of the template together with the constraint on the transformation introduce the shape prior. The deformation is guided by the maximization of the likelihood of foreground and background intensity models, so that we can refer to this approach as Competitive Deformation. In each region, the intensity is modelled as a smooth approximation of the original image. We represent the transformation using a Partition of Unity Finite Element Method, which consists in representing each component with polynomial approximations within local patches. A conformity constraint between the patches provides a way to control the globality of the deformation. We show several results on synthetic images, as well as on medical data from different modalities.

Introduction

Image segmentation is a fundamental topic in computer vision, which has motivated many works to cope with challenging issues such as noise, occlusions and low contrasted regions. A common approach is the introduction of prior knowledge in order to constrain the solution to remain close to a given class of shapes. Statistical models have been proposed, using for example Principal Component Analysis (PCA), like the well-known active shape models [START_REF] Cootes | Active shape models: Their training and application[END_REF]. Such techniques require careful training in order to capture the variability of the shapes, and do not enable segmentation if no training database is available.

Shape priors have also been incorporated in the level set framework via an additive shape term in the energy, that penalizes the dissimilarity between the level set function for segmentation and the one embedding the prior shape [START_REF] Paragios | Matching distance functions: A shape-toarea variational approach for global-to-local registration[END_REF][START_REF] Cremers | Towards recognition-based variational segmentation using shape priors and dynamic labeling[END_REF][START_REF] Chan | Level set based shape prior segmentation[END_REF][START_REF] Raviv | Prior-based segmentation by projective registration and level sets[END_REF]. In the model of Leventon et al. [START_REF] Leventon | Statistical shape influence in geodesic active contours[END_REF] and further works along the same line [START_REF] Tsai | Model-based curve evolution technique for image segmentation[END_REF][START_REF] Tsai | A shape-based approach to the segmentation of medical imagery using level sets[END_REF][START_REF] Bresson | A variational model for object segmentation using boundary information and shape prior driven by the mumfordshah functional[END_REF], a PCA of training shapes embedded in level set functions is computed in order to define a linear statistical model for the shape term. Non-linear versions have also been explored [START_REF] Cremers | Shape statistics in kernel space for variational image segmentation[END_REF]. For some applications, those models may however suffer from the uncontrollable topology changes allowed by the level set representation, which is most often not desirable when one wants to impose a shape prior.

A possible alternative to control the topology is to directly apply a geometric deformation to the prior shape, provided that the transformation is diffeomorphic [START_REF] Hong | Shape representation based on integral kernels: Application to image matching and segmentation[END_REF][START_REF] Saddi | Region-based segmentation via non-rigid template matching[END_REF]. Two-phase segmentation is performed by deforming a binary template towards the image. The prior is the template itself, and the shape constraint is conveyed through the choice of a class of allowed deformations. This choice is crucial: excessive constraint may result in poor segmentation results, whereas insufficient constraint may lead to a final shape too far from the prior.

The deformation can be guided by the maximization of the likelihood of userdefined intensity models, given the pixel values observed in both the foreground and the background. The approach bears similarities with Region Competition techniques [START_REF] Zhu | Region competition: Unifying snakes, region growing, and bayes/mdl for multiband image segmentation[END_REF], with the significant difference that the unknown variable is not the partitioning itself but the deformation of a predefined template. We will refer to this method as Competitive Deformation. Obviously, the choice of appropriate region intensity models is also essential. On one hand, piecewise-constant and global models are simple but their applicability is limited. On the other hand, piecewise-smooth and local models are more relevant in many cases, at the cost of an increased computational complexity.

We propose a variational formulation of two-phase, piecewise-smooth segmentation based on Competitive Deformation and a Partition of Unity Finite Element Method (PUFEM). The key idea is to represent the deformation field with polynomial approximations within overlapping local patches. The method includes a conformity constraint between the patches, which offers a good control over the range of the smoothness, and hence over the strength of the prior. Moreover, this framework provides efficient smooth representations of the region intensity models that naturally extrapolate beyond the boundary.

This paper is organised as follows. In section 2, we first give the basic formulation for two-phase segmentation by template competitive deformation. Then in section 3, we give a brief overview of the PUFEM framework before presenting our formulation of piecewise-smooth segmentation in section 4. In section 5, we show experimental results on illustrative synthetic images, and on medical images.

Basic formulation of Competitive Deformation

Let the open bounded set Ω ⊂ R d be the domain of a real-valued image I. Twophase segmentation aims at partitioning Ω into a foreground and a background that are homogeneous in terms of intensity properties. In the classical level set version of Region Competition [START_REF] Zhu | Region competition: Unifying snakes, region growing, and bayes/mdl for multiband image segmentation[END_REF], the optimal partition of Ω is obtained by evolving a level set function that embeds the region boundary, allowing undesirable topology changes in the case of prior-based segmentation. To control the topology of the result, we rather deform the characteristic function χ of a prior foreground region Σ:

χ(x) = 1 if x ∈ Σ 0 otherwise (1)
with a diffeomorphic geometric transformation ψ : Ω → ψ (Ω). The basic form of the two-phase competitive deformation problem reads:

min ψ,α1,α2    x∈Ω χ • ψ(x)r 1 (α 1 , x) + x∈Ω (1 -χ • ψ(x)) r 2 (α 2 , x) + γR (ψ)    (2)
where the r i : Ω → R are the region model functions that encode the intensity properties in the foreground and the background. Each r i usually depends on a set of parameters α i , e.g. the mean intensity value, the standard deviation, etc. R is the regularization constraint on ψ and γ a positive constant that controls its weighting. The minimization is carried out iteratively by alternating the following two steps:

(A) Considering the transformation ψ fixed, optimize and update the model's parameters α 1 and α 2 , (B) Considering α 1 and α 2 fixed, minimize w.r.t ψ.

Region intensity models

In this formulation, simple global statistics can be used to represent the region intensity properties, such as the well-known piecewise-constant case [START_REF] Chan | Active contours without edges[END_REF][START_REF] An | Region based image segmentation using a modified mumfordshah algorithm[END_REF], assuming Gaussian intensity distributions with known variance, i.e.:

r i (m i , x) = (I(x) -m i ) 2 (3)
m i being the mean intensity value in region i. The Gaussian assumption has practical limitations and is not valid in images showing more complex intensity distributions. In a similar template deformation context, Saddi et al. [START_REF] Saddi | Region-based segmentation via non-rigid template matching[END_REF] used the more general non-parametric model:

r i (p i , x) = -log p i (I(x)) (4)
where p i is the intensity probability density function in region i. Estimating global probability densities for the foreground and background using the whole image still has limitations in practice. Especially critical are the cases of cluttered and heterogeneous backgrounds and the presence of low-frequency artifacts such as illumination changes. Moreover, it may be difficult to obtain a precise positioning of the boundary since the local contributions of nearby pixels from both sides are diluted in the global estimation of the densities.

To overcome these limitations, it is natural to turn to local models, assuming smoothly-varying intensity distributions. In the Gaussian case with spacedependent mean values, this leads to piecewise-smooth segmentation. The error function reads:

r i (I i , x) = (I(x) -I i (x)) 2 + µ i |∇I i (x)| 2 (5)
where I i becomes a function that approximates I and is constrained to be smooth inside region i; µ i is the smoothing parameter. Vese and Chan [START_REF] Vese | A multiphase level set framework for image segmentation using the mumford and shah model[END_REF], and simultaneously Tsai et al. [START_REF] Tsai | Curve evolution implementation of the Mumford-Shah functional for image segmentation, denoising, interpolation, and magnification[END_REF], introduced a level-set formulation based on diffusion and on the work of Mumford and Shah [START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF]. This diffusion scheme requires the iterative resolution of a PDE with new boundary conditions at each update step (B), which is computationnally costly. Methods based on Gaussian convolution were recently proposed [START_REF] Mory | Fuzzy region competition: A convex two-phase segmentation framework[END_REF][START_REF] Brox | On the statistical interpretation of the piecewise smooth mumford-shah functional[END_REF][START_REF] Li | Implicit active contours driven by local binary fitting energy[END_REF], giving qualitatively similar results in a more efficient way. They enable the extrapolation of the model beyond the region boundary, to an extent that depends on the kernel's scale. In section 4, we present an alternative method based on finite elements with polynomial extrapolary properties.

Deformation models

Compliance with the shape prior is determined by the class of deformations allowed by the regularization constraint R. Therefore it has to be chosen carefully.

In particular, R must take into account the fact that a shape is invariant to some geometric, global transformations such as translation, rotation, scaling, shearing, etc. For example, several instances of a same shape are shown on Fig. 1. Despite Fig. 1. Several instances of a same shape their different positions, orientations and sizes, they all represent the same shape. Consequently, a relevant regularizer R shall not penalize such transformations. Within the template matching context, a similarity transformation is used in [START_REF] An | Region based image segmentation using a modified mumfordshah algorithm[END_REF] to segment synthetic images, excluding more complex deformations; in the non-rigid case, Saddi et al. [START_REF] Saddi | Region-based segmentation via non-rigid template matching[END_REF] constrain the deformation with a diffeomorphic fluid model, thus enabling the result to strongly deviate from the prior shape.

Methods based on local basis expansions of the deformation such as B-splines [START_REF] Kybic | Multidimensional elastic registration of images using splines[END_REF][START_REF] Rueckert | Nonrigid registration using free-form deformations: Application to breast MR images[END_REF] or Radial Basis Functions [START_REF] Bookstein | Principal warps: Thin-plate splines and the decomposition of deformations[END_REF][START_REF] Fornefett | Elastic registration of medical images using radial basis functions with compact support[END_REF] provide a suitable compromise between global and fluid models. As will be seen in the subsequent sections, we use a finite element registration framework based on a partition of unity, which enables to easily control the globality of ψ. Moreover, unlike Free Form Deformation methods, our regularization term does not penalize globally polynomial transformations.

Partition of Unity Finite Element Representation

In this section, we give an overview of the mathematical framework of the PUFEM [START_REF] Babuska | The partition of unity method[END_REF] that we use to represent a given scalar function: in our case, the approximation images I i in each region, and the components of the transformation ψ in each dimension. The basic idea is to locally fit the said scalar field with d-dimensional polynomials and smoothly blend them afterwards to obtain a regular representation.

Let F be a real-valued function defined on Ω. We define a set N of nodes distributed over Ω. A node n ∈ N is characterized by:

-a point c (n) ∈ Ω, called center of the node, -an open bounded subdomain Ω (n) ⊂ R d containing c (n) , called patch, -a function ϕ (n) : R d → R, called PU-function, -a set of ρ (n) functions B (n) = {p (n) r : Ω → R | r ≤ ρ (n) }, called the local basis.
We allow the patches to overlap and assume the families (Ω (n) ) n∈N and (ϕ (n) ) n∈N to fullfil the Partition of Unity conditions:

Ω ⊂ n∈N Ω (n) and ∀x ∈ Ω n∈N ϕ (n) (x) = 1 (6)
For the sake of computational efficiency, our nodes are distributed over a regular, rectangular array and each patch Ω (n) is a cuboid centered on c (n) . This configuration is illustrated on Fig. 2.a. The PU-function ϕ (n) has a compact support included in Ω (n) ; it is non-negative, equal to 1 at c (n) and vanishes with the distance to c (n) (cf. Fig. 2.b). Thus, any point x ∈ Ω belongs to 4 patches in 2D, and 8 patches in 3D, with weights given by the corresponding ϕ (n) (x). The basis functions p

(n) r
are the monomials of all degrees up to a user-defined maximum degree -e.g. in 2D and degrees up to 2: 1, x, y, x 2 , xy, y 2 -centered on c (n) , so that F is locally modelled at node n by a polynomial F (n) :

F (n) = r ρ (n) a (n) r p (n) r (7)
where the a (n) r are real coefficients. The global representation is then constructed by blending the F (n) with the PU-functions:

F = n∈N ϕ (n) F (n) = n∈N r ρ (n) ϕ (n) a (n) r p (n) r (8)
According to [START_REF] Tsai | A shape-based approach to the segmentation of medical imagery using level sets[END_REF], F (n) is as regular as the PU-functions per se. However, we want to impose a controllable, "long range" regularization, or rather, globality. To this end, we introduce the notion of non-conformity between two neighbouring nodes m and n through the energy:

S (m,n) κ ( F ) = Ω (m,n) |β| κ ϕ (m) ϕ (n) D β F (m) -D β F (n) 2 (9)
where β = (β 1 , β 2 , . . . , β k ) and D β is the partial derivative operator in the standard multi-index notations. This local energy has an intuitive interpretation: it penalizes F if its local representations at nodes m and n and their derivatives up to order κ differ in the overlapping region Ω (m,n) . The total conformity energy is then defined by:

S κ ( F ) = 1 2 n∈N m∈V(n) S (m,n) κ ( F ) (10)
where V (n) is the set of neighbours of node n in 4-connexity. This inter-node conformity constraint is a key feature of our method. It enables smooth representations of the region intensity models, that are naturally extrapolated beyond the boundary. As for the representation of the deformation, this energy is zero when all the local representations are equal, i.e. when ψ is globally polynomial. Thus, in the case of local affine bases, global translation, rotation, scaling and shearing are not penalized.

Piecewise-smooth segmentation with Competitive Deformation

We now detail our two-phase, piecewise-smooth segmentation formulation casting the Competitive Deformation problem into the PUFEM framework.

Our formulation

We address the deformation of the template characteristic function χ through its corresponding displacement vector field u = ψ -id, the components of which are represented as in ( 8) by sets of coefficients that we pile up into one vector a.

The regularization term R(ψ) in ( 2) takes the form of an inter-node conformity constraint S κ (a), defined as the sum of the constraints on each component of u.

In region i, I is approximated by a locally polynomial image I αi i , represented over a set of nodes N i by a set of coefficients α i according to ( 8) and ( 7):

I αi i = n∈Ni ϕ (n) i I (n) i
where

I (n) i = r ρ (n) i α (n) ir p (n) ir (11)
Most existing piecewise-smooth methods include regularization inside each phase (e.g. see ( 5)). Consequently, I i is not explicitly defined outside region i. We apply an inter-node conformity constraint S κi (α i ) involving all the nodes of N i , so that our energy functional reads:

E(a, α 1 , α 2 ) = Ω χ • ψ a (I -I α1 1 ) 2 + Ω (1 -χ • ψ a ) (I -I α2 2 ) 2 + γ S κ (a) + µ 1 S κ1 (α 1 ) + µ 2 S κ2 (α 2 ) ( 12
)
We minimize ( 12) by alternating steps (A) and (B) (see section 2) until convergence.

Step (A): Estimating the image approximations

Minimization w.r.t. the parameter set α i is achieved, considering the energy:

E i (α i ) = Ω χ i (I -I αi i ) 2 + µ i S κi (α i ) (13)
where we define χ 1 = χ • ψ and χ 2 = 1 -χ • ψ. The first term is a masked least square term that enforces I i to fit the image within region i. The second term is an inter-node conformity constraint that compels I i to be regular everywhere in Ω. In other words, I i results from a regularized approximation of I inside region i, and since no fitting constraint is imposed outside, it is extrapolated beyond the region border by regularization only. S κi is a quadratic function of the parameters α i and hence so is E i . Minimization is then achieved by classical linear regression. Saying that the derivatives of

E i w.r.t. the α (n)
ir vanish provides a system of linear equations which is sparse, since the nodes are only related to each other in 4-connexity, due to the overlap pattern between the patches (see Fig. 2.a). More precisely, we need to solve:

M i • α i = g i (14)
where M i is a symmetric, non-negative definite matrix of size n ρ

i , and g i a vector of length n ρ (n) i . Their entries are given in the appendix. We use a Conjugate Gradient descent, well-suited for solving sparse linear systems [START_REF] Shewchuk | An introduction to the conjugate gradient method without the agonizing pain[END_REF].

Step (B): Template Registration

We show that step (B) is equivalent to a classical registration problem based on a Sum of Square Difference criterion. Let χ be the signed characteristic function of the prior foreground Σ:

χ(x) = 1 if x ∈ Σ -1 otherwise (15)
Then we can replace χ by (1 + χ)/2 in ( 12). The region parameters α i being fixed, the energy to minimize w.r.t. a is:

- 1 2 Ω ( χ • ψ a ) r + γ S κ (a) (16)
where r = (I -

I α2 2 ) 2 -(I -I α1 1 )
2 . By writing:

χ • ψ a r = - 1 2 ( χ • ψ a -r) 2 -( χ • ψ a ) 2 -r 2 (17)
and since ( χ • ψ a ) 2 = 1 and r 2 is independent of the parameters a, we can reduce the energy to minimize to:

E u (a) = 1 4 Ω ( χ • ψ a -r) 2 + γ S κ (a) (18)
Therefore step (B) boils down to an SSD-based registration problem, with r being the reference and χ the template. Minimization of E u follows a global-tolocal strategy: we define a coarse-to-fine dyadic pyramid of node distributions.

Once the solution has been computed at one level, it is projected on the basis of the next finer level to provide an initialization. A detailed description of the minimization scheme is out of the scope of this paper and can be found in [START_REF] Makram-Ebeid | Non-rigid image registration using a hierarchical partition of unity finite element method[END_REF].

Results and discussion

We first present results on illustrative synthetic, noisy images (see Fig. 3). The first image, "χ", consists of two heterogeneous phases. The results show the robustness of our method to noise and strong intensity variations. The extrapolatory property is illustrated by figure (d) of the first line. The second image, "Treble Clef", consists of a homogeneous foreground on a heterogeneous background and contains occluding objects. The final segmentations are robust to occlusions and leaks.

We then apply our method to segment heart chambers on cardiac ultrasound (Fig. 4) and cine MR (Fig. 5) images. The final deformations are well-constrained so that the topologies of the prior shapes are preserved. On the MR images, the foreground region shall include the papillary muscles when segmenting the blood pool, which is challenging as they appear darker an may be confused with the myocardium.

"χ" and its prior shape "Treble Clef" and its prior shape On Fig. 6 we compare three deformation models: affine, PUFEM, and fluid. An affine deformation (a) is obviously too restricted to obtain an accurate segmentation from the prior shape that is used (Fig. 5.a bottom). On the opposite, the fluid model (c) allows strong local deformations and hence strong deviations from the prior shape. This is visible on the lower part of the right ventricle segmentation where we can see a curvature inversion when compared to the prior, and a consequent exclusion of the papillary muscles. The PUFEM model (b) provides a good compromise to have both accuracy and compliance with the prior shape.

Conclusion

We introduced a novel variational approach for two-phase, piecewise-smooth image segmentation based on prior shape Competitive Deformation. We cast our formulation into the Partition of Unity Finite Element framework, motivated by its long-range regularization properties, suitable for the class of transformations that is needed to abide by the shape prior. Indeed, the inter-node conformity constraint provides a good control over the globality of the deformation and does not penalize basic global transformations. This framework is also well-adapted for representing the region intensity models as smooth approximations of the original image on its whole domain. Our algorithm was successfully applied to challenging synthetic images and to medical images, with robustness to noise, occlusions and leaks. Appendix: entries of M i and g i

Estimating the image approximation in region i (step (A) in the minimization of the functional [START_REF] Saddi | Region-based segmentation via non-rigid template matching[END_REF]) boils down to the matricial equation:

M i • α i = g i (19)
where M i is a symmetric, non-negative definite matrix of size n ρ (n) i , and g i a vector of length n ρ (n) i . Their entries are given by:

m i (r,n)(s,m) = Ω ϕ (m) i ϕ (n) i   χ i p (m) is p (n) ir -µ i β κi D β p (m) is D β p (n) ir   + δ mn µ i Ω ϕ (n) β κi D β p (n) is D β p (n) ir (20)
g i (r,n) = Ω χ i ϕ (n) i p (n) ir I (21)
δ mn being the Kronecker delta equal to 1 if m = n and 0 otherwise. Since the support of a PU-function is included in the corresponding patch, m i (r,n)(s,m) equals zero when Ω (m) ∩ Ω (n) = ∅, i.e. when m and n are not neighbours in 4-connexity. Hence the sparseness of M i .

Fig. 2 .

 2 Fig. 2. Example of Partition of Unity configuration in 2D.

Fig. 3 .Fig. 4 .Fig. 5 .

 345 Fig. 3. Synthetic images "χ" and "Treble Clef": (a) Initializations. (b) Final segmentation. (c) Final piecewise-smooth approximations χI1 + (1 -χ)I2. (d) Final background approximations I2. (e) Final deformations.

Fig. 6 .

 6 Fig. 6. Comparison between three deformation models. Final segmentations are displayed on the first line and final deformations on the second line. (a) Affine model. (b) PUFEM model. (c) Affine + fluid model [12].


M Lopez-Martinez

C Canudas-De-Wit

F R Rubio

Stability of Asynchronous Feedback-Interconnected Dissipative Systems

This paper presents a method to analyze the stability of the feedback interconnection of a class of dissipative linear systems when the sampling associated to the feedback interconnection is asynchronous. We consider systems that are either Input/Output Strictly passive (IOSP), or systems which have bounded L2-gains less than one. The analysis is performed by using the concept of MAximum Sampling time preserving Dissipation (MASD), for each interconnected system. We investigate the impact of using the scattering transformation in the computation of the MASD, and we provide a numerical algorithm (based on a set of LMI's) that allows to choose the most suitable configuration for the interconnection.

I. INTRODUCTION

This paper presents a method to analyze the stability of the feedback interconnection of a class of dissipative linear systems when the sampling associated to the feedback interconnection is asynchronous. The considered scenario is shown in Fig. 1. These types of interconnections are present in many application fields such as remotely-operated systems [START_REF] Hokayem | Bilateral teleoperation: An historical survey[END_REF], interconnected vehicle control sub-systems, and more generally in component-based control design where synchronous exchange of information is not feasible. It is well known that, the continuous time interconnection of two systems in feedback form is stable when each system posses certain passivity property, or satisfies the small gain condition (see [START_REF] Van Der Schaft | L 2 -Gain and Passivity Techniques in Nonlinear Control[END_REF], [START_REF] Khalil | Nonlinear control systems[END_REF], [START_REF] Brogliato | Dissipative Systems Analysis and Control[END_REF]). As these two properties are fundamentally different (the first is phase shift sensible, while the second depends on its gain characteristics), it is expected that these conditions may behave differently when the output signals are sampled as shown in Fig. 1. In the framework of this work, the maximum sampling time preserving the stability of its original interconnection, may be different to each case.

Scattering transformation is a well known method whose change of variables is used to improve and/or to preserve the M. Lopez-Martinez and F.R. Rubio are with the Department of systems and automation, University of Seville, Spain mlm,rubio@esi.us.es C. Canudas-de-Wit is with the CNRS, at Gipsa-Lab, Grenoble, France carlos.canudas-de-wit@gipsa-lab.inpg.fr system original stability properties when the interconnecting signals are subject to delays. In addition, this transformation has been shown, in the continuous-time case [START_REF] Matiakis | The Scattering Transformation for Networked Control Systems[END_REF], [START_REF] Matiakis | Independent-of-Delay Stability of Nonlinear Networked Control Systems by Scattering Transformation[END_REF], [START_REF] Chopra | Delay-Independent Stability for Interconnected Nonlinear Systems with Finite L 2 Gain[END_REF], to imply a certain duality between passivity (IOSP) and L 2gain (with γ < 1) property. In these works, [START_REF] Matiakis | The Scattering Transformation for Networked Control Systems[END_REF], [START_REF] Matiakis | Independent-of-Delay Stability of Nonlinear Networked Control Systems by Scattering Transformation[END_REF], [START_REF] Chopra | Delay-Independent Stability for Interconnected Nonlinear Systems with Finite L 2 Gain[END_REF], the scattering transformation has been used in the context of continuous time system with fix or variable delays induced by the network, but the potential impact of an asynchronous sampling in the scenario shown in Fig. 2 has been scarcely studied. One of the objectives of this paper is to study how the scattering transformation impacts the maximum sampling time preserving the stability of its original interconnection. The lost of the original system properties such as dissipation and small gain, under the sampling process has been studied before. In [START_REF] Lin | Passivity and absolute stabilization of a class of discrete-time nonlinear systems[END_REF], [START_REF] De La Sen | Relationships between positive realness of continuos transfer functions and their digital counterparts[END_REF], [START_REF] De La Sen | Preserving positive realness through discretization[END_REF], [START_REF] Zhou | New characterization of positive realness and control of a class of uncertain polytopic discrete-time systems[END_REF] it has been shown that positive realness can be lost in the process of discretization, and that in order to preserve the property of passivity in sampled systems, it is necessary that the outputs of these systems depends on their inputs, forcing the system relative degree to be zero. In [START_REF] Hassibi | Control of asynchronous dynamical systems with rate constraints on events[END_REF] the authors have introduced a Lyapunov-based theory for asynchronous linear systems and shows how Lyapunov functions and controllers can be constructed by solving linear matrix inequalities (LMI). This theory was extended to nonlinear systems in [START_REF] De La Peña | Stability of nonlinear asynchronous systems[END_REF]. In [START_REF] Cueli | Passivity of Interconnected Asynchronous Discrete-Time Systems[END_REF], the interconnection of passive asynchronous systems was studied, and control design method (based in a set of LMI's) was introduced. The LMI's were constructed using an Euler approximation of the derivative of the Lyapunov function, being the approach different from [START_REF] Laila | Open and closed loop dissipation inequalities under sampling and controller emulation[END_REF], [START_REF] Nesic | Stabilization of sampled-data nonlinear systems via backstepping on their Euler approximate model[END_REF], [START_REF] Rohrer | Passivity considerations in stability studies of numerical integration algorithms[END_REF], where the approximation is introduced in the model.

This paper extends the results in [START_REF] Cueli | Passivity of Interconnected Asynchronous Discrete-Time Systems[END_REF], dealing only with passive system, to a more general dissipative systems, including L 2 -gain stable systems. The objective is also to devise a general method to find an optimal way to interconnect a dissipative system that tolerates the maximum possible sampled time intervals. The method includes the possibility to use or not the scattering transformation, and allows us to devise a numerical algorithm able to compute the Maximum sampling time preserving stability of such an interconnection.

The paper is organized as follows. In section II, we present the problem formulation for asynchronous feedbackinterconnected dissipative systems, and we also introduce some background on scattering theory with application to the interconnection of dissipative systems. The main results are given in Section III and IV. In Section III, we characterize the conditions required to preserve dissipation for asynchronous sampled systems, and in section IV, we provide necessary conditions to guarantee the stability of the interconnection. Our results are summarize in Section V.

II. PROBLEM FORMULATION AND PRELIMINARIES

Consider the feedback interconnection of two linear systems Σ 1 and Σ 2 , as shown in Fig. 1, with each system Σ i , of the following form:

Σ i : ẋi = A i x i + B i u i (1) y i = C i x i + D i u i . (2)
subject to a general dissipative condition given by

Vi (x, u) ≤ -y ⊤ i Q i y i + u ⊤ i R i u i + 2y ⊤ i S i u i -ρ i V i (x) (3)
where V i (x) is a positive definite function, that can be chosen for linear systems as V i (x) = x ⊤ i P i x i . We consider here two specific classes of systems characterized by the properties of matrices Q i , R i and S i (see [START_REF] Chellaboina | Exponentially Dissipative Nonlinear Dynamical Systems: A Nonlinear Extension of Strict Positive Realness[END_REF], [START_REF] Navarro-Lopez | Several dissipativity and passivity implications in the linear discrete time-setting[END_REF]): The problem considered in this paper is to find, for each system class C 1 , and C 2 the maximum sampling times, such that their dissipative and finite small gain original properties are preserved, and thus to ensure that the interconnection is stable. To this aim, we also consider the the possibility to use the scattering transformation in the setup shown in Fig. 2, and compare this possibility to the one without the transformation. This combination leads to four cases; A, B, C and D shown in Table I. There is a certain duality between these cases; a system that is IOSP (class C 1 ) has its L 2 -gain small than one after using the scattering transformation, S, and viceversa, a system with L 2 -gain small than one (class C 2 ) becomes IOSP after using the transformation S. This is discussed in detail next.

A. Passivity and L 2 -gain duality throughout scattering transformation

A way to analyze the stability of the interconnection of two systems is via Passivity or L 2 -gain. If both systems are passive or have their respective L 2 -gain equal or less than one, then the interconnection is stable. However, the two properties are fundamentally different; a system can preserve its passivity property while is L 2 -gain varies, while another can preserve its small gain property while the passivity is lost. This is illustrated by Fig. 3. Nevertheless, if the scattering transformation is used to interconnect the two sub-systems as in Fig. 2, i.e.

u = 1 √ 2 (U + J ), y = 1 √ 2 (U -J ) (4)
where u and y are respectively the input and the output of each system, and U and J are the ones resulting from this transformation, it has been shown (see Proposition 2.3.4 in [START_REF] Van Der Schaft | L 2 -Gain and Passivity Techniques in Nonlinear Control[END_REF]) that: if and only if a system is IOSP from u to y, then the resulting transformed system has a L 2 -gain less than one from U to J . This result can be easily understood following the next reasoning which will be useful for our next development.

Assume that we want the transformed system, Σ, to have L 2 -gain less than one. Therefore, it should hold that

V + ρV < 1 2 (U ⊤ U -J ⊤ J ) (5)
with ρ ≥ 0. After substitution of (4) in ( 5), we get

V + ρV -uy < 0.
This inequality is equivalent to say that there exists,for the system Σ, a positive definite function, w D (u, y),

w D (u, y) = y ⊤ Qy + u ⊤ Ru + 2y ⊤ Su > 0 (6)
where Q > 0 and R > 0, such that

V + ρV -uy ≤ -w D (u, y). (7)
To prove this, we substitute y = Cx + Du and define ξ T = [x T , u T ] in the above inequality, where

-w D (u, y) = -ξ ⊤ M ξ ≤ -εξ ⊤ ξ (8)
with

M = C ⊤ QC C ⊤ QD + C ⊤ S D ⊤ QC + S ⊤ C R + D ⊤ QD + D ⊤ S + S ⊤ D > 0,
and

ε = λ min (M ) > 0 ⇔ Q > 0, R > 0.
Finally, It is important to remark that always it is posible to find Q > 0 and R > 0 such that ε > 0 is small enough, which means that Σ must be an Input and Output Strictly Passive (IOSP) system.

B. Symmetry and modularity design in the scattering transformation

As there are many possibilities to define the scattering transformation, we wish to select one that make the transformation of one side independent to the other (modular design). In addition, and from the original problem, as defined in Fig. 1, the original system inter-connection requires that y 1 = u 2 and u 1 = -y 2 . Then, any possible transformation should then preserve this relation.

A good choice fulfilling such requirements is to define

U 1 = 1 √ 2 (u 1 + y 1 ), J 1 = 1 √ 2 (u 1 -y 1 ). (9)
for system Σ 1 , and

U 2 = - 1 √ 2 (u 2 + y 2 ), J 2 = - 1 √ 2 (u 2 -y 2 ) (10)
for system Σ 2 . It can be easily checked that this transformation give rise to the same transfer function,

Σ i : G i (s) = J i U i = u i -y i u i + y i = 1 -G i (s) 1 + G i (s) , (11)
in any of the two possible system coordinates. Note that the transformed system Σ i in terms of the new input U and the new output J can be written as

Σ i : ẋ = Āx + BU J = Cx + DU (12)
where

Ā . = A + B C, B . = B 1 √ 2 I + D , C . = - √ 2(I + D) -1 C and D . = (I + D) -1 (I -D).
The dissipative properties the transform system can be expressed, using (3), as

V ≤ -J ⊤ QJ + U ⊤ RU + 2J ⊤ SU -ρV, (13)
where the new dissipative matrixes are defined as

Q . = 1 2 (Q + R -2S) R . = 1 2 (Q + R + 2S) S . = 1 2 R ⊤ -Q ⊤ . (14)
III. CHARACTERIZATION OF THE MAXIMUM SAMPLING TIME PRESERVING DISSIPATION: MASD In this section we aim to characterize the maximum sampling time that preserve the system original dissipation properties. We will treat the general case of the dissipation property including the two class of systems C 1 , and C 2 , with and without the scattering transformation described in the previous section. The results will be presented as the solution of a LMIs guaranteing that the sampled system will satisfy the corresponding dissipative condition.

A. Discrete-time system description

We have to re-formulate our problem in the discrete time framework starting with the discrete time representation of the system Σ i in (1),

Σ i,k : x k+1 = A k x k + B k u k (15) y k = Cx k + Du k , (16)
where A k = e AT k and B k = A -1 (A k -I)B. and the system Σ i in [START_REF] Hassibi | Control of asynchronous dynamical systems with rate constraints on events[END_REF],

Σ i,k : x k+1 = Āk x k + Bk U k J k = Cx k + DU k . (17)
Next, the equivalent definition of the dissipative property needs to be reformulated also in this discrete-time framework. One possible way is to make a discrete-time approximation of (3), or [START_REF] De La Peña | Stability of nonlinear asynchronous systems[END_REF] according to the considered case, by using the approximation

V ≈ V k+1 -V k T k ,
where V k = x ⊤ k P x k , and proceeding as in previous section to write the associated discrete-time LMI, which can be used to find the maximum value for T k , preserving its dissipative property.

Note that each of the two considered classes of system can be then evaluated by the proper substitution of the associated (Q, R, S, ρ) values.

B. Original system coordinates

The dissipation inequality (3) evaluated in the setup shown in Fig. 1, gives the following LMI

A ⊤ k P T k A k -P T k A ⊤ k P T k B k B ⊤ k P T k A k B ⊤ k P T k B k + + C ⊤ QC + ρP C ⊤ (QD + S) (D ⊤ Q + S ⊤ )C R + D ⊤ QD + D ⊤ S + S ⊤ D ≤ 0. (18)
which in a compact notation writes as W (T k ) + Z ≤ 0, with W (T k ) being the first matrix sampling dependent, and Z the matrix capturing the particular dissipative property of the original system. Note that here the values of (Q, R, S, ρ) , and P are given from the original system properties, and then we only need to search for the maximum value of T k that verifies the LMI. Note also that from continuity of the solution, (there exist always a solution for infinitely small T k ), a limit solution for T * k > 0 will always exist. The definition of MASD, T * k > 0, follows then from this observation, i.e.

T * k = max T k >0 {T k : W (T k ) + Z ≤ 0} (19)

C. Using the Scattering transformation

Our second option will be to perform a similar analysis in the transformed system coordinates after applying the previously presented scattering transformation, in the setup given by Fig. 2. Proceeding as before, but now using (13), we get

Ā⊤ k P T k Āk -P T k Ā⊤ k P T k Bk B⊤ k P T k Āk B⊤ k P T k Bk + + C⊤ Q C + ρP C⊤ ( Q D + S) ( D⊤ Q + S⊤ ) C R + D⊤ Q D + D⊤ S + S⊤ D ≤ 0. (20)
which in a compact notation writes as W (T k ) + Z ≤ 0. The definition of MASD T * k > 0, follows as before, i.e.

T * k = max T k >0 T k : W (T k ) + Z ≤ 0 (21)
Remark 1: Note that both formulations ( 18) and (20) yield similar forms of the LMIs, but with matrices that inside have different values. Therefore, the expected values for each of the MASD; T * k , and T * k will be different.

D. Examples

1) First-order IOSP system: This example compares the MASD obtained for a dissipative system Σ of class C 1 (Case A in Tab. I), with and without the scattering transformation. For this, consider the following input-output strictly passive (IOSP) system,

G(s) = 1 s + 1 + d
and recall that the use of the scattering transformation, makes the transformed system Σ have L 2 -gain less than one (Case D in table I). Figure 4 shows the evolution of the T * k and the T * k for different values of the parameter d. It can be noticed that the resulting MASD is larger when the scattering transformation is used ( T * k > T * k ). The example reveals how the use of the scattering transformation improves over the one without transformation. 2) Underdamped Second Order System: This example compares the MASD obtained for a dissipative system Σ of class C 2 (Case B in Tab. I), with and without the scattering transformation. Consider the system

G(s) = Kω 2 n s 2 + 2δω n s + ω 2 n ( 22
)
with parameters K = 0.865, ω n = 6 and δ = 0.5. This system Σ is not passive, but has L 2 -gain equal to 0.999. The second option is to use the transform system using the scattering matrix. Because the duality indicated in previous section, the transform system Σ belongs to class C 2 (Case C in Tab. I). Figure 5-(b) shows the λ = min(Re{G(jω)}) as a function of the sampling time. It can be observed that λ > 0 for all sampling time, which shows that the best option (in terms of stability) to interconnect this system is by the scattering transformation. One of the paper goal is to select the formulation that gives the larger MASD, while preserving the stability of the interconnection. Stability issues under such interconnection framework are studied next.

IV. STABILITY PROPERTIES OF ASYNCHRONOUS SAMPLED SYSTEMS

Let us assume that each system1 is sampled using a zero order hold, and as before, that the sampling time-intervals T 2,k , and T 1,k are not constant but are multiple-integers. For simplicity we consider the case in which

T 2,k = n k T 1,k , with n k ∈ {1, 2, . . . , n max }, but similar results hold if n k is of the form n k ∈ {1, 1/2, . . . , 1/n max }.
A. Stability conditions in the original system coordinates Proposition 1: Consider two dissipative systems2 Σ i , fulfilling each one its respective dissipative inequality [START_REF] Chellaboina | Exponentially Dissipative Nonlinear Dynamical Systems: A Nonlinear Extension of Strict Positive Realness[END_REF], and interconnected asynchronously as shown in Fig. 1. Assume that there exists for each system, Σ i , a MASD T *

i , i = 1, 2, verifying T 2,k = n k T 1,k ≤ min{T * 1 , T * 2 }
, Then, the following holds:

• the feedback asynchronous interconnection of two dissipative systems is asymptotically stable if the following LMI is satisfied

Q 1 + R 2 S 1 -S ⊤ 2 S ⊤ 1 -S 2 R 1 + Q 2 > 0, (23)
• If ρ i > 0, i = 1, 2, then the feedback interconnection is also exponentially stable, • For the two classes C 1 , and C 2 of considered systems, the condition (23) holds, and hence the interconnection is asymptotically stable.

Proof: Assuming that there exist

T * i , i = 1, 2, verifying T 2,k = n k T 1,k ≤ min{T * 1 , T * 2 }
, such that each sub-system Σ i fulfills its respective dissipative inequality, implies that that both systems must verify a dissipative property for the smaller sampling time, T 1,k , that is

Σ 1 : V 1,k+1 -V 1,k ≤ ≤ -T 1,k (y ⊤ 1,k Q 1 y 1,k + u ⊤ 1,k R 1 u 1,k + 2y ⊤ 1,k S 1 u 1,k + ρ 1 V 1,k ) Σ 2 : V 2,k+1 -V 2,k ≤ ≤ -T 1,k (y ⊤ 2,k Q 2 y 2,k + u ⊤ 2,k R 2 u 2,k + 2y ⊤ 2,k S 2 u 2,k + ρ 2 V 2,k )
As both systems are feedback interconnected it follows that u 2,k = y 1,k and y 2,k = -u 1,k . Taking as a Lyapunov function candidate,

V k = V 1,k + V 2,k , the rate of change of this function ∆V k = V k+1 -V k is given as ∆V k = V 1,k+1 -V 1,k + V 2,k+1 -V 2,k ≤ (24) -T 1,k (y ⊤ 1,k Q 1 y 1,k + u ⊤ 1,k R 1 u 1,k + 2y ⊤ 1,k S 1 u 1,k + ρ 1 V 1,k + +u ⊤ 1,k Q 2 u 1,k + y ⊤ 1,k R 2 y 1,k -2u ⊤ 1,k S 2 y 1,k + ρ 2 V 2,k ) < 0 introducing ρ . =min(ρ 1 , ρ 2 ) yields, ∆V k ≤ -T 1,k y ⊤ 1,k (Q 1 + R 2 )y 1,k + u ⊤ 1,k (R 1 + Q 2 )u 1,k + 2y ⊤ 1,k (S 1 -S ⊤
2 )u 1,k + ρV k < 0 It follows straightforward that if ρ = 0 asymptotic stability can be reached if the two systems are ZSD, and the condition (23) hold, insuring that ∆V k ≤ 0.

The second item in the proposition is simple to check by observing that if ρ > 0, then ∆V k is an strictly decreasing function, and hence the interconnection become exponentially stable.

Finally, the last item is easy to check by simple substitution of the definition of the (Q i , R i , S i , ρ i ) in the equation (23).

The proposition can be particularized for every kind of dissipative systems, i.e. in terms of passivity or in terms of L 2 -gain. The proposition is also valid for systems interconnected with the scattering transformation, as shown next.

B. Stability conditions using scattering transformation

Proposition 2: Consider two dissipative systems 3 Σi , fulfilling each one its respective dissipative inequality (20), and interconnected asynchronously by means of a scattering transformation as shown in Fig. 2. Assume that there exists for each system, Σi , a MASD T * i , i = 1, 2, verifying T2,k = n k T1,k ≤ min{ T * 1 , T * 2 }, Then, the following holds:

• the feedback asynchronous interconnection of two dissipative systems is asymptotically stable if the following LMI is satisfied

Q1 + R2 S1 -S⊤ 2 S⊤ 1 -S2 R1 + Q2 > 0, (25)
• If ρ i > 0, i = 1, 2, then the feedback interconnection is also exponentially stable, • For the two classes C 1 , and C 2 of considered systems, the condition (25) holds, and hence the interconnection is asymptotically stable.

Proof: The proof is similar to the one in Proposition 1, and then it is omitted.

Remark 2: Substituting ( 14) in (25) the LMI can be expressed in terms of the original dissipative properties of the systems with no scattering transformation.

Remark 3: The condition given by the two previous propositions, can be also used in an inverse manner to find for the most appropriate maximum sampling time ensuring stability. Besides, it may be possible also to use this analysis for control design along the same lines than in the works of [START_REF] Cueli | Passivity of Interconnected Asynchronous Discrete-Time Systems[END_REF].

C. Example

Consider the interconnections of two systems of the form (22) as described in the Example III-D.2. The goal of this example is to evaluate by simulation the interconnection behavior for the two possible scenarios shown in Fig. 1 and Fig. 2, with the addition of a reference step signal applied to the system Σ 1 to better visualize the system transient responses.

The first simulation tests are done by the proper selection of the sampling time validating the stability conditions: 3 Both systems are assumed to be zero state detectable (ZSD). The second example, shown in Figure 7, repeat the previous simulation but now with T k = 0.55 s. This sample time selection implies that (T * k < T k < T * k ). As the computation in the Example III-D.2 have shown, the use of the scattering matrices allows for larger sampling interval for stability, but also an improvement in the system transient behavior.

T k =0.3 s, which implies that T k < T * k < T * k .

V. CONCLUSIONS

In this paper we have presented a method to analyze the stability of the feedback interconnection of a class of linear systems when the sampling is asynchronous. The analysis has been made taking into account dissipatedness properties of the original continuous-time system.

We have characterized specific conditions in terms of LMIs allowing to compute the larger sampling time interval preserving the original system dissipative properties. For that, we have also considered the possibility of using the scattering transformation. This transformation allows, in some cases, to improve the system stability over the situation where it is not used, because longer sampling time intervals can thus be tolerated.

We have presented stability conditions for the asynchronous interconnection of two general dissipative systems. The conclusions are extended to the case of using the scattering transformation in the interconnection. Finally, we have introduced a numeric procedure for choosing an optimal way to interconnect a dissipative system.

Fig. 1 .

 1 Fig. 1. Asynchronous Feedback Interconnection. System Σ 1 and system Σ 2 are feedback interconnected throughout the asynchronous, and possible time-varying, sampling T 1,k , and T 2,k .

Fig. 2 .

 2 Fig.2. Interconnection by means of scattering. Scattering may be used as a mean to improve the system stability under asynchronous sampled interconnection.

C 1 : 2 :

 12 Input/Output Strictly Passive systems (IOSP);Q = δI > 0, S = -1 2 I and R = ǫI > 0 C Finite small L 2 -gain systems (γ < 1); Q = I, S = 0 and R = -γ 2 I

Fig. 3 .

 3 Fig. 3. Output passive system with different L 2 -gains (left). System with different passivity properties preserving its L 2 -gain (right).

Fig. 4 .

 4 Fig. 4. Comparison of MASD for a first-order IOSP system with and without scattering transformation

Figure 5 -

 5 (a) shows the L 2 -gain versus the sampling time. It can be observed that the gain over-exceeds one at T * = 0.46.

Fig. 5 .

 5 Fig. 5. Evolution of dissipative properties. (a)-Evolution of the L 2 -gain γ of G(s) as a function of the sampling time, and (b)-Evolution of λ = min(Re{G(jω)}) as a function of the sampling time.

Figure 6 -

 6 (a) shows the system outputs corresponding to the interconnection in Figure. 1 (Case B in Tab. I), and Figure 6-(b) the system output using the scenario shown in Figure 2 (Case C in Tab. I), including the scattering transformation. It can be seen how the inclusion of the scattering transformation increases the achieved performance.

Fig. 6 .

 6 Fig.6. Time-evolution of the system output y 1 (t). Sampling time fulfilling the stability conditions for both cases: (a) without the use of the scattering transformation, and (b) using the scattering transformation.

Fig. 7 .

 7 Fig.7. Time-evolution of the system output y 1 (t) with T * k < T k < T * k . The case of sampling time not fulfilling the stability conditions is shown in (a) which does not the use of the scattering transformation. The case (b) using the scattering transformation is stable and has a good performance.

TABLE I DIFFERENT

 I CASES TO BE STUDIED IN THE PAPER.

		System class	System class
		C 1	C 2
	without scattering	A	B
	with scattering	C	D

The notion of system here designs, either Σ i , or Σi , according with the considered case.

For completeness we also assume that both systems are zero state detectable (ZSD).

ACKNOWLEDGMENTS

The authors gratefully acknowledge MEC (Ministerio de Educación y Ciencia español) for founding this work under grants PR2007-0282 and DPI2007-64697, and CNRS (Centre National de la Recherche Scientifique, France) and INRIA-RA for allowing that the stay of the first author was possible.


Maxim Fradkin

Cybèle Ciofolo

Benoit Mory

Gilion Hautvast

Marcel Breeuwer

Comprehensive Segmentation of Cine Cardiac MR Images

A typical Cardiac Magnetic Resonance (CMR) examination includes acquisition of a sequence of short-axis (SA) and long-axis (LA) images covering the cardiac cycle. Quantitative analysis of the heart function requires segmentation of the left ventricle (LV) SA images, while segmented LA views allow more accurate estimation of the basal slice and can be used for slice registration. Since manual segmentation of CMR images is very tedious and time-consuming, its automation is highly required. In this paper, we propose a fully automatic 2D method for segmenting LV consecutively in LA and SA images. The approach was validated on 35 patients giving mean segmentation error smaller than one pixel, both for LA and SA, and accurate LV volume measurements.

Introduction

A typical Cardiac Magnetic Resonance (CMR) examination includes acquisition of a sequence of short-axis (SA) and long-axis (LA) images covering the cardiac cycle. Quantitative analysis of the heart function (ejection fraction, myocardial thickness and thickening, etc.) requires segmentation of the left ventricle (LV) SA images. LA views, in addition, allow more reliable LV shape estimation, especially near the apex and the valve plane. The later can be used for the accurate and reproducible determination of the basal SA slice, which is known to be one of the major factors of interobserver variability in LV measurements [START_REF] Weaver | Magnitude and causes of interobserver discrepancies in CMR volumme measurements: critical importance of choice of the basal slice[END_REF][START_REF] Anderson | Normal cardiac magnetic resonance measurements and interobserver discrepancies in volumes and mass using the papillary muscle inclusion method[END_REF]. Automatic multi-view segmentation methods, enabling to obtain highly reproducible LV measurements [START_REF] Van Geuns | Automatic quantitative left ventricular analysis of cine MR images by using threedimensional information for contour detection[END_REF], are thus highly desirable.

Many publications propose automatic or semi-automatic methods for segmenting the LV in SA CMR images [START_REF] Mitchell | Multistage hybrid active appearance model matching: segmentation of left and right ventricles in cardiac MR images[END_REF][START_REF] Kaus | Automated segmentation of the left ventricle in cardiac MRI[END_REF][START_REF] Niessen | Geodesic deformable models for medical image analysis[END_REF][START_REF] Paragios | A variational approach for the segmentation of the left ventricle in cardiac image analysis[END_REF][START_REF] Jolly | Automatic segmentation of the left ventricle in cardiac MR and CT images[END_REF] or in multiple views [START_REF] Lelieveldt | Multi-view active appearance models for consistent segmentation of multiple standard views: application to long-and short-axis cardiac MR images[END_REF][START_REF] Lotjonen | Statistical shape model of atria, ventricles and epicardium from short-and long-axis MR images[END_REF][START_REF] Van Assen | SPASM: A 3D-ASM for segmentation of sparse and arbitrarily oriented cardiac MRI data[END_REF], to mention a few. Since SA CMR images, acquired over multiple breath-holds, are often misregistered due to patient motion or inconsistent respiration, 3D segmentation and analysis methods require a registration preprocessing step [START_REF] Moore El | A high resolution dynamic heart model based on averaged MRI data[END_REF][START_REF] Lotjonen | Correction of movement artifacts from 4-d cardiac short-and long-axis MR data[END_REF][START_REF] Perperidis | Statistical shape model of atria, ventricles and epicardium from short-and long-axis MR images[END_REF][START_REF] Hautvast | Alignment of short axis and long axis cine cardiac MR images[END_REF][START_REF] Slomka | Automated registration of multiple single breath-hold cardiac MRI images[END_REF]. Alternatively to image-based alignment, one can use for this purpose LV contours extracted from multiple views [START_REF] Li | Registration of short and long-axis images in cine cardiac MRI[END_REF]. In contrast to [START_REF] Van Geuns | Automatic quantitative left ventricular analysis of cine MR images by using threedimensional information for contour detection[END_REF], where LA contours had to be drawn manually to initialize SA segmentation, we propose a fully automatic method following similar workflow.

The paper is organized as follows: Section 2 gives a concise description of the proposed method; in Section 3 validation results are presented and discussed; conclusions are given in Section 4.

Among 2D segmentation methods, explicit active contours [START_REF] Kass | Snakes: active contour models[END_REF] remain popular, since the underlying discrete contour representation is sufficient to capture the fairly smooth shape of the LV. However, this representation often produces slowly converging algorithms that are quite sensitive to initial conditions and parameterization. While there are various solutions to tackle those issues within the original framework, an interesting alternative representation using B-splines has been proposed in [START_REF] Menet | BSnakes: implementation and application to stereo[END_REF]. This continuous representation is compact and allows the semi-local control of smooth curves. Those advantages have motivated the design of our spline-based deformable template of the LV.

Most active contour methods are based on the minimization of an objective criterion and tend to be locally trapped by spurious image features during the optimization process. The greedy optimization framework [START_REF] Williams | A fast algorithm for active contours and curvature estimation[END_REF] used in this paper offers a valuable compromise between computational complexity, robustness and flexibility. Not only it gives in most cases an acceptable solution by avoiding local minima, but also allows easy incorporation of dissimilar image constraints and prior knowledge (e.g., both differentiable and non-differentiable). We follow a two-stage approach: First, we robustly detect LV position on LA view(s) using standard image acquisition geometry and rigid template deformation. The precise delineation of the LV is done by deforming the template locally, through greedy optimization. Second, the segmented LA contours are used for initialization of deformable templates in the relevant SA slices, which are then adjusted using the same, as for LA, optimization paradigm.

Deformable template geometry

Our parametric template relies on interpolating splines controlled by as few nodes as possible. More specifically, the myocardium is modeled as a closed ribbon structure, composed of an imaginary centerline (dashed in Fig. 1) and a variable width. Both the centerline C(s) = [x(s), y(s)] and the ribbon width w(s) are continuous spline interpolations of a discrete set of {p k = (x k , y k , w k )} samples defined at each node. Among the advantages of this compact representation is the natural coupling between the endocardium and the epicardium (the inside and outside contours). We also define two enclosed regions M and B, corresponding respectively to the myocardium and blood pool areas, as shown in Fig. 1.

Minimization criterion

The optimal solution of our segmentation problem is the set of template parameters minimizing a criterion that expresses the goodness of fit between the geometrical model and the image evidence, given some a priori knowledge. In the remainder, I is the image and the λ i 's are scalar weights balancing the various terms. Let p = p k = (x k , y k , w k ) T , k ∈ 1, N be our parametric model. The problem can now be formalized as follows:

min p F (p, I) = F s (C, w) shape + F c (C i , C o , I) contour + F r (M, B, I)
region Shape: Shape constraint is composed of three terms:

F s (C) = λ 0 F regularity (C, w) + λ 1 F smoothness (C) + λ 2 F similarity (C) (1)
The first one regularizes template width:

F regularity (C, w) = 1 0 |w ′ (s)|ds
The smoothness term controls centerline curvature:

F smoothness (C) = 1 0 |κ(s) -κ 0 |ds
where κ(s) stands for the curvature of the centerline C(s) and κ 0 is a desired curvature ,e.g. setting it to average curvature constrains the contour to be circular (see Sec. 2.5).

The similarity term constraints the template curve shape to be affinely similar with a pre-defined shape C. It is defined as an error of the best affine transformation T between the given contour C and a pre-defined contour C:

F similarity (C) = 1 0 |C(s) -T ( C(s))| 2 ds
Contour: The contour term stands for the contour contrast, both local (first two terms) and global (last term):

F c (C i , C o , I) = λ 3 1 0 ∇I(C i (s)).n(s)ds-λ 4 1 0 |∇I(C o (s)).n(s)| ds+λ 5 (m-b) (2)
where ∇I is the image gradient, n(s) contour's normal, m and b are the average intensities of the regions M and B, respectively.

Region: The myocardium and the blood pool gray levels, except for the papillary muscles region, should be homogeneously distributed. Therefore we can define F r using a minimal variance criterion such as:

F r (M, B, I) = λ 6 |M| M |I(x, y) -m| 2 dxdy + λ 7 |B| B I(x, y) -b 2 dxdy (3)
where region M (resp. B) has average intensity m (resp. b) and area |M| (resp. |B|). In order to increase robustness to the expected presence of papillary muscles, which breaks the blood pool homogeneity assumption, a rather effective approach is to remove the contribution of a certain percentage (e.g., 15%) of the darkest pixels in region B from the variance computation.

Optimization strategy

We use a greedy optimization scheme embedded in a coarse-to-fine approach to optimize simultaneously the nodes position and ribbon width. Due to using very few nodes and a global criterion, each optimization step influences an important portion of the model, thus allowing potential jumps over local minima.

Unfortunately, casting normals with such a sparse sampling makes the solution dependent on the initial position of the nodes, hence decreasing parameterization invariance. We compensate for this drawback of the original scheme by a systematic sliding of all nodes along the centerline between each iteration.

Segmentation of long axis images

Template creation and initial pose: The initial template should resemble the LV shape, so its centerline can be approximated by an ellipse which axes are aligned with the main axes of the LV; it is slightly bended inside for modeling valve shape (Fig. 2-left). For initialization of LA template, we need to define its position (p c ) and orientation ( -→ u ). These parameters can be first obtained in patient coordinates, and then easily transformed into the image ones. When two LA views are available, we reasonably assume (due to standard acquisition planning) that an intersection of the two LA image planes is approximately aligned with the real LA of the heart (Fig. 2-center) and therefore can be taken as -→ u (it is taken co-oriented with the Z axis of the SA). The template origin p c is obtained as an intersection of the axis -→ u with the middle slice of the SA volume. When only one LA view is available, we align -→ u with Z axis of the SA volume, and the template origin is defined as a closest -→ u point to the center of SA volume p o , since a typical SA acquisition is centered on LV (Fig. 2-right).

Template deformation: The deformation is done in two steps: coarse (rigid) and fine (local). At the coarse step, the initial template is affinely optimized using only two minimization terms, namely contour (Eq. 2) and region (Eq. 3), since the shape is implicitly conserved by the transformation. Note also that in the contour term only global contrast contributes at this stage (λ 3 , λ 4 →0). This step, allowing to robustly detect the LV position, is followed by a local deformation, using all three terms. In the contour constraint, only local contrast is now put to contribution (λ 5 →0), while in the shape constraint (Eq. 1) only regularity and similarity terms are used (λ 1 →0). In addition, contour resulting from the coarse step is used as C, which preserves certain contour similarity with the initial shape and in particular allows to softly constrain valve concavity.

Segmentation of short axis images

Template creation and initial pose: The SA template centerline can be approximated by a circle. Thus, SA template creation is quite obvious since one only needs to define the template center and radius (Fig. 3-left). Given the LA segmentation results, the LV LA contour(s) intersections are found at each SA slice. If two LA contours are available, the template parameters are found by fitting a best circle through 4 intersection points (Fig. 3-center), while if only one contour is given, the template center is defined as the middle of the intersection segment and its radius as a half of its length (Fig. 3-right). Note that if no LA images are available, the initialization of SA templates can still be done using robust ring detection algorithm [START_REF] Breeuwer | Automatic segmentation of short-axis cine cardiac magnetic resonance[END_REF].

Template deformation: The SA template deformation follows the same coarseto-fine scheme as the LA. The optimization terms used at coarse step are exactly the same as in Sec. 2.4. At the fine step, the contour and the region terms are again the same as for LA, while the shape constraint comprises only regularity and smoothness terms (in the later one, κ 0 is set to the average contour curvature to impose contour circularity, thus also substituting the similarity term).

Validation

Data: The clinical data used for validating our method included 35 patients. SA cine CMR scans comprised 9-14 slices and 15-50 phases; all SA images were 256x256 in size, covering a field of view ranging from 360x360 mm up to 450x450 mm. For 5 patients two LA views (2 chambers (CH) and 4CH) were available, while for the remaining ones only 4CH LA views were acquired. All LA images were 256x256 in size, covering a field of view ranging from 350x350 mm up to 460x460 mm. Manual delineations, done by an expert, were available for all enddiastolic (ED) images (325 in total of SA and 40 for LA).

Results: To compare the result of our automatic segmentation method with the manually delineated contours, we measured the mean distance (ε mean ), the rootmean-square (RMS) distance (ε rms ) and the maximum distance (ε max ) between the resulting contours and the manual ones, using method explained in [START_REF] Hautvast | Automatic contour propagation in cine cardiac magnetic resonance images[END_REF]. For SA images, mean positioning errors was 1.3 mm for the endocardium and 1.5 mm for the epicardium, for a pixel size ranging from 1.4 mm to 1.8 mm (see Table 3). For LA images, mean positioning errors was 1.3 mm for the endocardium and 1.1 mm for the epicardium, for a pixel size ranging from 1.3 mm to 1.8 mm. An excellent correlation was found between "manual" and "automatic" volumes (r =0.90, P < 0.001), with mean difference -10 ± 18 ml (P = 0.002). Note that all results were obtained with exactly the same method settings and without any user interaction. An example of segmentation results for one patient (2 LA views and 8 SA slices) is shown in Fig. 4. In our current implementation, the overall segmentation time is around 1 second per image on a 3.2 GHz PC, thus segmentation of a typical image phase (SA+LA) takes about 10 seconds. Discussion: LA template initialization proved to be robust even using only one LA image, the same stands for the initialization of SA templates from only one LA segmentation contour. Note also that SA initialization was robust to slices misregistration, since possible initialization inaccuracies were compensated by the rigid deformation. Producing a mean positioning error inferior to one pixel for both SA and LA images, comparable with inter-expert drawing variability [START_REF] Van Assen | SPASM: A 3D-ASM for segmentation of sparse and arbitrarily oriented cardiac MRI data[END_REF], the method proved accurate and robust. We validated the method for ED phase only, since the segmentation of the whole cardiac cycle can be then obtained using automatic contour propagation [START_REF] Hautvast | Automatic contour propagation in cine cardiac magnetic resonance images[END_REF], shown to preserve its accuracy within acceptable ranges while being much faster than phase-by-phase segmentation.

Conclusion

We developed a fully automatic, robust and accurate method for segmentation of CMR SA and LA acquisitions, enabling easy introduction of prior knowledge and constraints in the optimization framework. The method was extensively validated by comparing the segmentation results with manually drawn endocardial and epicardial contours for 35 patients (325 SA slices and 40 LA images). The positioning error was inferior to one pixel for both SA and LA images and is comparable with interobserver variability. The resulting contours can be used for SA slice alignment, contour propagation and accurate LV parameter estimation.

Fig. 1 .

 1 Fig. 1. Template geometry (left); Detailed zoom (center).

Fig. 2 .

 2 Fig. 2. LA template: centerline as an ellipse with upper bending modeling the valve (left); LA template initial pose using SA volume, represented as a cube, and two LA views, represented as planes (center); Similarly to later, LA template initial pose with SA volume and only one LA view (right).

Fig. 3 .

 3 Fig.3. SA template centerline as a circle (left); SA template initial pose using intersection of two LA contours with a SA volume slice (center); SA template initial pose using intersection of only one LA contour with a SA volume slice (right).

  (mm) εrms (mm) εmax (mm) εmean (mm) εrms (mm) εmax (mm) SA 1.27 ± 0.54 1.56 ± 0.70 3.29 ± 1.62 1.56 ± 0.85 1.93 ± 1.04 4.06 ± 2.18 LA 1.32 ± 0.41 1.69 ± 0.55 4.55 ± 1.84 1.13 ± 0.42 1.52 ± 0.67 4.01 ± 2.27Table1. Endocardial and epicardial positioning errors compared to manual contours.

Fig. 4 .

 4 Fig. 4. LV segmentation example for one patient (from left to right, up to bottom): LA 2CH image; LA 4CH image; consecutive slices of SA volume from apex to valve.


Behavior of a poroelastic material bonded onto a vibrating plate is investigated in the low frequency range. From the analysis of dissipation mechanisms, a model accounting for damping added by the porous layer on the plate is derived. This analysis is based on a 3-D finite element formulation including poroelastic elements based on Biot displacement theory. First, dissipated powers related to thermal, viscous and viscoelastic dissipation are explicited. Then a generic configuration (simply supported aluminium plate with a bonded porous layer and mechanical excitation) is studied. Thermal dissipation is found negligible. Viscous dissipation can be optimized as a function of air flow resistivity. It can be the major phenomenon within soft materials, but for most foams viscoelastic dissipation is dominant. Consequently an equivalent plate model is proposed. It includes shear in the porous layer and only viscoelasticity of the skeleton. Excellent agreement is found with the full numerical model.

INTRODUCTION

Porous materials like polymer foam and glass wool are widely used for noise control in several engineering activities such as aeronautics and automotive industries. Their properties are twofold : sound absorption and damping of the nearby structure.

Excited through the coupling with the fluid [1] or directly by the structure when in contact with it [2], the skeleton participate to damping through viscoelastic dissipation.

In the literature, the efficiency of porous material is shown in applications dealing with sound absorption, transmission loss of panels [3,4] or noise in enclosure coupled with elastic panels [START_REF] Atalla | Aeroacoustics of flight vehicles : theory and practice[END_REF]6]. Some are based on finite element calculations including poroelastic elements [3,[START_REF] Atalla | Aeroacoustics of flight vehicles : theory and practice[END_REF]. In these studies damping effect is shown for various boundary conditions and materials, but no general tendencies is drawn about the nature of dissipation. Other reports give analytical predictions [START_REF] Biot | [END_REF]8] of the importance of viscous dissipation induced by the vibration of the skeleton, but they are approximative and limited to 1-D or 2-D applications.

The present analysis is based on a partition of dissipated and reactive powers in poroelastic media, as similarly proposed by Rasolofosaon [9]. In order to treat any 3-D structure, the calculation relies on a 3-D finite element formulation [10] with poroelastic elements based on Biot displacement theory [11].

Well suited for low frequency range analysis, this formulation has been validated for various applications [10,[START_REF] Dauchez | Experimental validation of 3-D poroelastic finite elements[END_REF]. Note that a similar study can be done using the mixed displacement-pressure (u, P ) formulation [13]. However, this latter formulation is not used here, since the presented work has been initiated with the (u, U ) displacement model. Both formulations have been shown to be equivalent [14]. However the cost of (u, U ) formulation is higher since three fluid displacement components are involved while the (u, P ) formulation relies on a single pressure degree of freedom for the fluid phase.

A comparison of the different terms of dissipated power gives the relative importance of dissipation mechanisms : viscous, thermal or viscoelastic. The analysis is performed on a structures comprising a simply supported plate associated with a porous layer excited by a point force (figure 1). The influence of its thickness, stiffness and air flow resistivity is studied. The goal is to determine major dissipation mechanisms in order to derive a simplified model of this two layers structrure. This model gives the characteristics of an equivalent plate accounting for mass, stiffness and viscoelastic damping added on the plate. Well suited for low frequency analysis, this model provides an efficient alternative to the full Biot model in terms of reduction of memory requirement and computation time.

POWER PARTITION

The discretized motion equation of the two layer system is given by Panneton

[10] [Z e ] 0 0 [Z p ] {w} u U = {F} , (1)
where [Z e ] and [Z p ] are respectively the impedance matrix related to the elastic and poroelastic media. {F } is the vector of nodal forces , applied to the elastic and poroelastic media. {w}, { u U } are the complex amplitude of the nodal displacements of the plate, the skeleton and the air comprised in the pores respectively. In the case where the porous material is bonded onto the plate, there is a continuity of the displacements at the interface

u = w, U n = u n , (2)
where subscript n denotes normal displacement. Tangential fluid displacement is not affected by continuity relations at the interface. These linear relations are applied by multiplication with a contraction matrix. Harmonic time dependence of the form e jωt is assumed so that instantaneous quantities are expressed as

a(t) = ℜ(a e jwt ), (3)
where a is complex amplitude and ℜ denotes real part. Instantaneous input power P(t) of the discretized system is the product of the instantaneous velocities by the instantaneous input forces F(t),

P(t) = ẇ(t), u(t), U (t) F(t) , (4)
where symbols , { } and ˙denote respectively line and column vector, and time derivation. This instantaneous power can be decomposed in two components related to power absorbed by the structure and energy exchanged between the excitation and the structure, such as

P(t) = D 1 + cos2(ωt) + R sin2(ωt). (5)
The mean power dissipated during one cycle D and the amplitude of the reactive power R derive from the complex power P , given by

P = 1 2 jω w * , u * , U * F (6) = D + jR, (7)
where * denotes a complex conjugate quantity. Using equation (1), complex power can be rewritten in the form

P = 1 2 jω w * , u * , U * [Z e ] 0 0 [Z p ] {w} u U , (8)
depending only on the uncoupled impedance matrix, and on the solution of the global system accounting for interface continuity. For the plate, impedance matrix

[Z e ] = [K e ] -ω 2 [M e ], (9)
comprises a real mass matrix [M e ] and a complex stiffness matrix [K e ] accounting for various loss phenomena, such as structural dissipation, acoustic radiation and loss through boundary conditions. Impedance matrix of the poroelastic media is given by Panneton

[10] [Z p ] = [K p ] + jω[C p ] -ω 2 [M p ]. (10)
This matrix comprises a real mass matrix [M p ], a complex viscous loss matrix [C p ] and a complex stiffness matrix

[K p ] = [K ss ] [K sf ] [K sf ] T [K f f ] , (11)
where [K ss ] is the complex stiffness matrix of the solid phase, accounting for structural dissipation. [K sf ] and [K f f ] are complex matrix related respectively to elastic coupling and bulk modulus of the fluid phase. They include also thermal losses from fluid to skeleton. Coupling, viscous and thermal losses are fonction of φ porosity, α ∞ tortuosity, σ air flow resistivity, Λ and Λ ′ viscous and thermal characteristic lengths respectively, according to Johnson-Allard-Lafarge theory [11, [START_REF] Johnson | Theory of dynamic permeability and tortuosity in fluid-saturated porous media 16[END_REF]16,17]. Relations between low and high frequency effects are governed by shape factors. Viscous shape factor M is defined by

M = 8µα ∞ σφΛ 2 , ( 12
)
where µ is the viscosity of the air. Thermal shape factor M ′ is defined by

M ′ = 8k ′ 0 φΛ 2 , ( 13
)
where k ′ 0 is the thermal permeability. This quantity is deduced from equation 13 assuming M ′ =1, as in a cylindrical pore.

Using equations (8-11), complex power can be split into several terms, such as

P = D ke + D ks + D kf + D cp + j(R ke + R kp + R me + R mp ), (14)
where R ke and R kp are elastic reactive powers of the plate and porous media. R me and R mp are inertial reactive powers of the plate and porous media. D ke , D ks , D kf and D cp are dissipated powers related to initial damping of the plate, structural, thermal and viscous dissipations within the porous material. Their expressions are

R ke = 1 2 ℑ jω w * K e {w} , (15)
R kp = 1 2 ℑ jω u * , U * K p u U , (16)
R me = 1 2 ℑ -jω 3 w * M e {w} , (17)
R mp = 1 2 ℑ -jω 3 u * , U * M p u U -jω 2 u * , U * C p u U , (18)
D ke = 1 2 ℜ jω w * K e {w} , (19)
D ks = 1 2 ℜ jω u * K ss {u} , (20)
D kf = 1 2 ℜ jω u * , U * 0 [K sf ] [K sf ] T [K f f ] u U , (21)
D cp = 1 2 ℜ -jω 2 u * , U * C p u U , (22)
where ℑ and ℜ denotes respectively imaginary and real parts.

A global loss factor η g can be calculated for the whole structure. It is expressed as the ratio of total dissipated power over the elastic reactive power of the structure

η g = D ke + (D ks + D kf + D cp ) R ke + R kp . ( 23
)
This expression is consistent with modal loss factor.

DISSIPATION MECHANISMS

Dissipation mechanisms are studied for a simply supported aluminium plate with a bonded academic porous layer. Continuity of displacement is insured between the plate and the porous layer according to equation (2).

Along the other faces of the porous layer, skeleton and fluid displacements are not constrained : this means that external fluid loading is neglected. This hypothesis, usually assumed for a heavy structure immerged in a light fluid, might have an influence on our conclusions [18-19]. However, plate-foam experimental tests [20] have shown good correlation between experiments and numerical simulations with the above assumption.

The rectangular plate (22 cm × 28 cm) is meshed by 13 × 13 thin shell quadragular elements. Excitation is achieved by a normal point force on the plate, located at (7.7 cm, 9.8 cm). The porous layer is meshed by 13 × 13 × 7 linear hexahedric elements. This mesh is suitable for the present study, but could be reconsidered for a more critical application due to slow convergence of linear poroelastic elements [START_REF] Dauchez | Rayonnement acoustique des structures -Vibroacoustique, Interactions fluide-structure[END_REF].

Calculations performed for various porous materials and thicknesses show that only structural and viscous dissipations are relevant. Thermal dissipation is always found less than 2 % of the total dissipation. This shows that the fluid is not significantly compressed : this is due firstly to the nature of excitation that does not create a straight compression of the fluid in the pore like an acoustical excitation would do, and secondly to the free boundary conditions of the porous layer that do not constrain the fluid in a given volume.

The relative importance of structural and viscous dissipation is now analyzed by varying the related parameters. Structural dissipation induced by the porous layer on the plate is function of the structural loss factor η 2 and the Young's modulus E 2 of the skeleton. It is studied for two couples (E 2 , η 2 ), including most of mineral wool or polymer foam characteristics : (E 2 =60 kPa, η 2 =0.07) used for a soft material and (E 2 =400 kPa, η 2 =0.15) for a stiff and more dissipative material. The aluminium plate properties are : thickness h 1 =1 mm, Young's modulus E 1 =69 GPa, Poisson ratio ν 1 =0.33 and density ρ 1 = 2778 kg/m 3 . The loss factor of the plate η 1 is set to 0.01 to account for various dissipation phenomena such as structural damping, acoustic radiation and losses through boundary.

In the low frequency range, viscous dissipation is mainly related to flow resistivity σ of the material. Its influence is determined from 10 3 to 10 8 Nsm -4 . Characteristic lengths varies with resistivity so that the used material remains realistic : viscous shape factor M is kept constant (Equation 12) and the ratio between thermal and viscous characteristic lengths is fixed to 2.5 . The other parameters are fixed and correspond to usual foam or glass wool characteristics : porosity φ=0.98, tortuosity α ∞ =1.3, skeleton density ρ 2 =40 kg/m 3 and Poisson ratio ν 2 = 0. The influence of thickness h 2 of the porous layer is also investigated, being set to 2 cm or 5 cm.

Figure 2.a presents viscous dissipation relative to total dissipation, D cp /(D ks + D kf + D cp ), as a function of flow resistivity, for the two thicknesses and materials. It is given for the first bending mode (65 Hz < f 11 < 68 Hz) of the structure.

An optimal value of flow resistivity appears for each case. Its location tends to smaller flow resistivity when thickness increases. This is consistent with the calculation of Okuno [8], given for a sealed poroelastic beam undergoing bending deformation

σ ≃ φ K 0 2.5 h 2 2 f i , (24)
where K 0 is the bulk modulus of the air and f the frequency. Relative viscous damping increases with the thickness. It reaches a maximum of 84% for the soft material, being the major dissipating mechanism. However it is limited to 27% for the stiff material : structural damping is the major dissipating mechanism. Figure 2.b presents global loss factor of the multilayer, η g . For very low air flow resistivity, this factor depends only on the structural damping and on the thickness of the porous layer. For a thickness of 2 cm, only stiff polymer foam adds a significant amount of damping. For soft material, damping becomes significant for a thicker layer and only where the viscous dissipation is important : a maximum loss factor of 6% is then reached with 5 cm. However, according to Okuno [8], the optimized value of air flow resistivity depends on the frequency : it is not suitable for a large frequency range.

The most reliable way to get damping in the low frequency range is to optimize viscoelastic dissipation with a most dissipative and stiff skeleton. In that case, the amount of added damping is no more dependent on air flow resistivity. Then, air flow resistivity can be then optimized for other purposes such as acoustic absorption.

EQUIVALENT PLATE MODEL

Viscoelastic dissipation is now assumed to be the major source of damping. In this context, only the skeleton behavior is relevant. The porous layer can be considered as a monophasic viscoelastic media and characteristics of an equivalent plate to the two layer structure can be determined. Development of such an equivalent plate is interesting to reduce computational cost. Moreover, this model can be coupled with an acoustical admittance model [22] when the porous layer is connected to a cavity : in such a case, both acoustical and structural effects of the porous layer, on the cavity and on the elastic structure, are taken into account respectively.

The geometry of th problem is described by figure 3. The two layers are characterized by Young's modulus E 1 and E 2 , Poisson's ratios ν 1 and ν 2 , densities ρ 1 and ρ 2 and thicknesses h 1 and h 2 . Subscript 1 refers to the plate and 2 to the porous layer.

For pure bending deformation of the structure, the equivalent plate parameters can be calculated from Ross-Kerwin-Ungar (RKU) theory [START_REF] Ross | Structural Damping[END_REF] of multilayer plate. Total bending rigidity D 12 is simply the sum of bending rigidities D 1 and D 2 of the two layers related to the neutral fiber of the plate,

D 12 = D 1 + D 2 , (25)
where D 1 and D 2 are given in the appendix by equations ( 48) and (49). The equivalent density is given by

ρ 12 = ρ 1 h 1 + ρ 2 h 2 h 1 + h 2 . ( 26
)
Equivalent loss factor is given by the ratio of the imaginary and the real part of the total bending rigidity D 12 . However and as expected RKU theory overestimates the equivalent loss factor in comparison with numerical calculation, for a relatively thick layer (figure 4.b). This confirms that porous layer, past a certain thickness, do not exhibit a pure bending deformation.

A more accurate model is achieved by taking into account shear strain in the porous layer (figure 3). Hamilton's principle is used to get the equation of motion of the two layer plate.

Classical thin plate theory is used for the plate. Only calculation related to the porous layer are exposed. For the porous layer, displacement vector {u} T = u, v, w of a particle is assumed to be

   u(x, y, z, t) = -z ∂w ∂x -(z -z 1 )ψ x , v(x, y, z, t) = -z ∂w ∂y -(z -z 1 )ψ y , w(x, y, t) = w, ( 27
)
where w is the deflection of the plate, ψ x and ψ y deviation angles due to shear strain (figure 3). The components of the strain tensor {ǫ 2 } T = ε x , ε y , γ xy , γ xz , γ yz are

             ε x = -z ∂ 2 w ∂x 2 -(z -z 1 ) ∂ψx ∂x , ε y = -z ∂ 2 w ∂y 2 -(z -z 1 ) ∂ψy ∂y , γ xy = -2z ∂ 2 w ∂x∂y -(z -z 1 )( ∂ψx ∂y + ∂ψy ∂x ), γ xz = -ψ x , γ yz = -ψ y . (28)
Stress-strain relations are given by

{σ 2 } = [H 2 ]{ǫ 2 } ( 29
)
where {σ 2 } = σ x , σ y , τ xy , τ xz , τ yz is the stress tensor and [H 2 ] is the viscoelastic tensor for a bidimensional isotropic material

[H 2 ] = E 2 1 -ν 2 2       1 ν 2 ν 2 1 1-ν 2 2 1-ν 2 2 1-ν 2 2       . ( 30
)
Strain energy V of the multilayer is the sum of bending strain energy of the plate, bending and shear strain energies of the porous layer

V = V bending 1 + V bending 2 + V shear 2 (31)
where,

V b1 = 1 2 E 1 1 -ν 2 1 S z 1 -z 1 ε x , ε y , γ xy   1 ν 1 ν 1 1 1-ν 1 2      ε x ε y γ xy    dxdydz, (32)
V b2 = 1 2 E 2 1 -ν 2 2 S z 2 z 1 ε x , ε y , γ xy   1 ν 2 ν 2 1 1-ν 2 2      ε x ε y γ xy    dxdydz, (33)
V s2 = 1 2 S z 2 z 1 γ xz , γ yz τ xz τ yz dx dy dz. (34) (35)
The kinetic energy of the system is given by

T = 1 2 S m 12 ẇ2 dS, (36)
with the equivalent density

m 12 = ρ 1 h 1 + ((1 -φ)ρ s + φρ 0 )h 2 , ( 37
)
where φ is the porosity, ρ s is the density of the solid comprising the skeleton and ρ 0 is the density of the air. Applying Lagrange's equations respectively to each variable, w, ψ x and ψ y , and summing the equations related to ψ x and ψ y gives the two equations of motion

D 1 ∆∆w + D 2 ∆∆w + D 4 ∆θ + m 12 ∂ 2 w ∂t 2 = 0, (38)
D 3 ∆θ + D 4 ∆∆w -C 2 θ = 0, (39)
with θ = ∂ψx ∂x + ∂ψy ∂y and ∆ = ∂ 2 ∂x 2 + ∂ 2 ∂y 2 , and where D 3 and D 4 are bending rigidities given in appendix by equations ( 50) and (51). For a steady state motion, w and θ are assumed to be of the form

w(x, y, t) = w 0 sin(k x x) sin(k y y) sin(ωt), (40) θ
(x, y, t) = θ 0 sin(k x x) sin(k y y) sin(ωt), (41)
where k x and k y are wavenumbers associated to the directions x and y. They are related to the wavenumber k by

k 2 = k 2 x + k 2 y . (42)
If the plate is simply supported and a × b sized, its modes (r a , r b ) correspond to

k x = π r a a and k y = π r b b , (43)
where r a and r b are modal orders along x and y directions respectively (r a > 1, r b > 1). The dispersion equation of the system (38 -39) can then be written in the form

ω 2 = k 4 D 12 m 12 C s (k 2 ), (44)
where C s (k 2 ) is the correction factor of the bending rigidity D 12 , accounting for shear, given by

C s (k 2 ) = 1 - D 2 4 k 2 D 12 (D 3 k 2 + C 2 ) . (45)
As C s is less than 1, shear will lower the resonance frequencies of the system. C s depends on the wavenumber, thus on frequency. For a simply supported plate, C s can be calculated for discrete values of frequency, corresponding to natural frequencies of the plate, given by equations ( 43) and (44). However C s can be expressed continuously as a function of frequency by solving equation ( 44) that is a polynomial equation of the third degree on k 2 .

The developed equivalent model and the full numerical model using poroelastic elements are compared in figure 4, in the case described in section 3 where viscoelastic dissipation is dominating. This is also the case where shear is the most important in the porous layer due to its large thickness. The porous layer is of a 5 cm thick and corresponds to the stiff material (E 2 =400 kPa, η 2 =0.15) with air flow resistivity σ= 5000 Nsm -4 . Characteristics of the equivalent plate are first calculated for the given configuration and frequency range. Then they are introduced as properties of a single plate. Equivalent loss factor is given by the ratio of the imaginary and the real part of the corrected bending rigidity C s (k 2 ) D 12 .

The RKU model is found insufficient for the plate-foam configuration : damping and eigenfrequencies are overestimated. The proposed equivalent model gives a good estimation of mean quadratic velocity of the plate and global loss factor for the first three modes. Comparison has not been performed for higher modes because of the computational cost of the poroelastic finite element method. However it can be seen that the equivalent plate model tends to underestimate damping as frequency increases. One reason is the increase of viscous forces along with frequency. This shows the limitation of the present model to the low frequency range.

CONCLUSION

Partition of dissipated and reactive powers is presented according to the 3-D formulation coupling elastic and poroelastic elements. Dissipation analysis has been performed in the low frequency range for a two layers structure comprising a porous layer bonded onto a plate. It shows that viscous and viscoelastic dissipation dominate. According to Okuno [8], viscous dissipation can be optimized by choosing a proper flow resistivity : it becomes the major dissipation mechanism within soft materials. However for stiff polymer foams, viscoelastic dissipation in the skeleton is widely dominating. Because viscous dissipation requires tuning of flow resistivity as a function of thickness and frequency, optimizing damping is most efficiently achieved by using the most viscoelastic and stiff porous material. Determination of an optimized flow resistivity is rather connected to sound absorption properties.

Consequently, an equivalent plate model accounting for the effect of the porous layer on the plate has been derived : it includes shear in the porous layer and only viscoelasticity of the skeleton. Good results have been obtained in comparison with the model using poroelastic elements for the three first modes, even for a thick porous layer.

The analysis performed on the two layers structure demonstrated the validity of the method based on power partition. Relying on a 3-D finite element formulation, this approach can be applied to any structure associating acoustic, elastic and poroelastic media. With a best understanding of the behavior of such structures, suitable simplified models can be derived avoiding the use of poroelastic elements, well-known for high computational cost.

APPENDIX

V b1 = D 1 2 S ∂ 2 w ∂x 2 2 + ∂ 2 w ∂y 2 2 + 2ν 1 ∂ 2 w ∂x 2 ∂ 2 w ∂y 2 + 2(1 -ν 1 ) ∂ 2 w ∂x∂y 2 dS, (46)
V b2 = 1 2 S D 2 ∂ 2 w ∂x 2 2 + ∂ 2 w ∂y 2 2 + 2ν 2 ∂ 2 w ∂x 2 ∂ 2 w ∂y 2 + 2(1 -ν 2 ) ∂ 2 w
with the bending rigidities (51)

D 1 = E 1 1 -ν 2 1 z1 -z1 z 2 dz = E 1 h 3 1 12(1 -ν 2 1 ) , (48)
D 2 = E 2 1 -ν 2 2 z2 z1 z 2 dz = E 2 1 -ν 2 2 z 3 2 -z 3 1 3 , (49)
D 3 = E 2 1 -ν 2 2 z2 z1 (z -z 1 ) 2 dz, (50)
D 4 = E 2 1 -ν 2
V s2 = 1 2 C 2 S ψ 2 x + ψ 2 y dS, (52)
with the shear rigidity

C 2 = κ h 2 E 2 2(1 + ν 2 ) , (53)
and κ, accounting for the variation of the shear stresses and strains through the thickness, taken to 5 6 .

  -z 1 ) dz.

Figure 1 .Figure 2 .

 12 Figure 1. Simply supported plate with a bonded porous layer and mechanical excitation. Figure 2. Percentage of viscous dissipation relative to total dissipation in the porous layer (a), global loss factor of the multilayer (b), as a function of flow resistivity for two materials : long dashed line --, 2 cm thick stiff material; continuous line -, 5 cm thick stiff material; dotted line • • • , 2 cm thick soft material; short dashed line ---, 5 cm thick soft material. Figure 3. Geometry of the two layer plate with shear.Figure 4. Comparison of equivalent model results with full discretized structure for a 5 cm thick low resistive (σ=5000 Nsm -4 ) and stiff (E 2 =400 kPa, η 2 =0.15) foam layer. (a) Mean quadratic velocity of the plate , (b) Global loss factor of the multilayer : continuous line -, plate with poroelastic elements; dashed line ---, equivalent plate with pure bending; dotted line • • • , equivalent plate with shear.

Figure 4 .

 4 Figure 1. Simply supported plate with a bonded porous layer and mechanical excitation. Figure 2. Percentage of viscous dissipation relative to total dissipation in the porous layer (a), global loss factor of the multilayer (b), as a function of flow resistivity for two materials : long dashed line --, 2 cm thick stiff material; continuous line -, 5 cm thick stiff material; dotted line • • • , 2 cm thick soft material; short dashed line ---, 5 cm thick soft material. Figure 3. Geometry of the two layer plate with shear.Figure 4. Comparison of equivalent model results with full discretized structure for a 5 cm thick low resistive (σ=5000 Nsm -4 ) and stiff (E 2 =400 kPa, η 2 =0.15) foam layer. (a) Mean quadratic velocity of the plate , (b) Global loss factor of the multilayer : continuous line -, plate with poroelastic elements; dashed line ---, equivalent plate with pure bending; dotted line • • • , equivalent plate with shear.

Figure 1

 1 Figure 1.

Figure 2

 2 Figure 2.

Figure 4 .

 4 Figure 4.

ACKNOWLEDGMENTS

The authors are grateful to PSA Peugeot Citroën and Agence de l'Environnement et de la Maîtrise de l'Energie (Ademe) who have supported this study.


Nicolas Dauchez

Sohbi Sahraoui

Noureddine Atalla

Convergence of poroelastic finite elements based on Biot displacement formulation

come