Romina Spalazzese
email: romina.spalazzese@di.univaq.it

Towards Mediator Connectors for Application level Interoperability

Keywords: D.2.12 [Software Engineering]: Interoperability Design, Theory Mediator connectors, protocol interoperability

My research work aims at overcoming the heterogeneity barriers that prevent networked systems to interoperate, thus enabling them to communicate responding to the continuously evolving networked environment. The thesis focusses on the concept of mediator connector as the key enabler for system interoperability.

The proposed work about mediator connectors is investigated in the scope of the Connect 1 European project.

INTRODUCTION

Our everyday activities are increasingly dependent upon the assistance of digital systems that pervade our living environment. Key technologies such as the Internet, the Web, and wireless computing devices and networks are now calm technologies in the sense of Marc Weiser's definition [START_REF] Weiser | The computer for the 21st century[END_REF]. They can indeed be qualified as ubiquitous, even if converged computing and networking technologies have still not reached the maturity envisioned by the ubiquitous computing and subsequent pervasive computing and ambient intelligence paradigms.

However, the current ubiquity of digital systems is technology-dependent. The efficacy of integrating and composing networked systems is proportional to the level of interoperability of the systems' respective behaviors, from application-down to network-layers. This leads to a landscape of islands of networked systems, although interoperability bridges may possibly be deployed among them. Further, the fast pace at which technology evolves at all abstraction layers increasingly challenges the lifetime of networked systems in the digital environment.

My research work aims at overcoming the heterogeneity barriers that prevent networked systems to interoperate, thus enabling them to communicate responding to the continuously evolving networked environment. The next section further discusses the interoperability issue and the concept of mediator connectors as a possible solution. Section 3 briefly introduces the core concepts of the approach and Section 4 concludes summarizing the related challenges.

THE INTEROPERABILITY PROBLEM

Middleware positions itself as the architectural paradigm enabling to effectively network together heterogeneous systems, specifically providing upper layer interoperability. Middleware bridges the gap between application programs and the lower-level hardware and software infrastructure in order to coordinate how application components are connected and how they interoperate, especially in the networked environment. As a matter of fact, middleware is yet another technological block, which also creates islands of networked systems.

Interoperable middleware have been introduced to overcome middleware heterogeneity. However, the solutions remain rather static, requiring either the use of a proprietary interface or a priori implementation of protocol translators. In general, interoperability solutions solve protocol mismatch among middleware at the syntactic level, which is too restrictive. This is even truer when one considers the many dimensions of heterogeneity, including software, hardware and networks, which now arise in ubiquitous networking environments, and that require fine tuning of the middleware according to the specific capacities embedded within the interacting parties. Thus, interoperable middleware can at best solve protocol mismatches arising among middleware aimed at a specific domain. Indeed, it is not possible to a priori design a universal middleware solution that will enable effective networking of digital systems, while spanning the many dimensions of heterogeneity now arising in networked environments and which will also increase dramatically in the future.

Mediator then stands as a core architectural paradigm for today's and future systems that increasingly need be connected. The mediator concept was early introduced to cope with the integration of heterogeneous data sources [START_REF] Wiederhold | Mediators in the architecture of future information systems[END_REF][START_REF] Wiederhold | The conceptual basis for mediation services[END_REF]. However, with the significant development of Web technologies and given abilities to communicate openly for networked systems, many heterogeneity dimensions shall now be mediated among which protocol mediation that is concerned with behavioral mismatches that may occur during interactions.

A key challenge for today's systems architectures is to embed the necessary support for automated mediation, i.e., the connector concept needs to evolve towards the one of mediator connectors. Indeed, the actual systems with which communications will take place cannot be anticipated at design time due to today's open networking and further continuous evolution of networked systems. As such, connectors not solely coordinate the interaction behaviors of connected systems but also mediate those behaviors to enable actual interactions.

Automated mediation of heterogeneous protocols basically relies on: (i) the adequate modeling of processes abstracting the behavior of the protocols to be bridged (ii) the definition of a matching relationship between the process models that sets the conditions under which protocol interoperability is supported (iii) the elicitation of an algorithm that computes an appropriate mapping between matching process models.

THE APPROACH

A number of solutions to automated protocol mediation have recently emerged. Although proposed algorithms manipulates formally grounded process models, most solutions are discussed informally, making it difficult to assess their respective advantages and drawbacks. They further remain rather vague on the definition of enforced matching relationship.

What is needed is a new and formal foundation for mediator connectors from which: protocol matching and associated mapping relationships may be rigorously defined and assessed. These relationships may be automatically reasoned upon, thus paving the way for on the fly synthesis of mediator connectors. To the best of my knowledge such an effort has not been addressed in the past.

The problem of protocols mediation concerns the interoperability between two protocols. To explain the problem, let us consider for example two instant messengers: Windows Messenger(W ) and Jabber Messenger(J). With respect to protocol meditation we are interested to formalize: i) how to identify whether the messenger protocols share similar intent ii) how to synthesise the mediator connector that enables them to interoperate despite of mismatches. In more detail with "protocols share similar intent" we mean that given the interaction protocols PW and PJ of W and J respectively, part of their behavior is complementary thus showing an interaction potentiality. Thus, we expect to find, at a given level of abstraction, similarities in the structure of the protocol representations. This leads to formally analyze such alike protocols to find, if it exists, a suitable mediator that allows the interoperability that otherwise would not be possible. I assume to know the behavioral specifications BW and BJ of the protocols PW and PJ respectively and also to have an ontology mapping between the alphabets of the mismatching protocols that defines a common alphabet. The high level process that I am investigating to solve this problem is made by the following steps. The first one deals with the transformation of BW and BJ into canonical forms, CBW and CBJ respectively, in order to simplify them and reason on minimized protocols. The second step deals with the check of the compatibility between CBW and CBJ . If they are compatible then one tries to synthesise the mediator connector to allow PW and PJ to interoperate otherwise the interoperability is not possible.

The compatibility check (also called matching), if succeeds, singles out the parts of the structures of the protocols that are expected to be similar, i.e., the parts of actual communication between PW and PJ . Then a mapping will build the behavioral representation for the abstract mediator for such part of protocols. The subsequent synthesis produces the actual implementation of the mediator.

To validate the proposed approach, being protocols rigorously defined, it can be formally proved the correctness and the completeness of the synthesised mediator.

CONCLUSION

With the networked environment being increasingly open and dynamic, the seamless composition and related connection of systems becomes a prime requirement. However, the openness of the environment comes along with great heterogeneity in the networked systems, which challenges their connection. Mediator connector then appears as a paradigm of choice to effectively overcome behavioral mismatches among the protocols run by the networked systems that need to coordinate. However, while the mediator paradigm has deserved much attention over the last few years, including research towards supported automated mediation, key principles remain loosely defined.

The proposed approach aims at formally characterizing the mediator connector problem. First results in this direction are presented in [START_REF] Spalazzese | A formalization of mediating connectors: Towards on the fly interoperability[END_REF].

A set of challenges arise in the area of software systems engineering (that would be interesting for future work), from theoretical foundations to specify the interaction behavior of networked systems to run-time methods and tools to turn specifications into running protocols, and vice versa.

ACKNOWLEDGMENTS

The work is partly supported by: the project Connect No 231167 of the Future and Emerging Technologies (FET) programme within the ICT theme of the Seventh Framework Programme for Research of the European Commission and the Italian PRIN D-ASAP. I would also like to thank Scuola Superiore Sant'Anna and ISTI CNR for presently hosting me after the Earthquake.


Samuel Kadoury

Nikos Paragios

Décembre

Inferring Preoperative Reconstructed Spine Models to Volumetric CT Data through High-Order MRFs

Keywords: Registration, physical modeling, image segmentation, articulated 3D spine model, Higher-order MRFs Recalage d'images, modélisation physique, segmentation, modèle 3D articulé de la colonne vertébrale, MRF haut-niveau

In this paper, we introduce a novel approach based on higher order energy functions which have the ability to encode global structural dependencies to infer articulated 3D spine models to CT volume data. A personalized geometrical model is reconstructed from biplanar X-rays before spinal surgery in order to create a spinal column representation which is modeled by a series of intervertebral transformations based on rotation and translation parameters. The shape transformation between the standing and lying poses is then achieved through a Markov Random Field optimization graph, where the unknown variables are the deformations applied to the intervertebral transformations. Singleton and pairwise potentials measure the support from the data and geometrical dependencies between neighboring vertebrae respectively, while higher order cliques are introduced to integrate consistency in regional curves. Optimization of model parameters in a multi-modal context is achieved using efficient linear programming and duality. A qualitative evaluation of the vertebra model alignment obtained from the proposed method gave promising results while the quantitative comparison to expert identification yields an accuracy of 1.8 ± 0.7mm based on the localization of surgical landmarks.

Introduction

Deformable anatomical models are powerful tools for recovering the shape of a patient's anatomy when only partial information or sparse image data is available. For example in orthopedic surgery, 3D computer-generated models have assisted specialists in surgical planning and instrument navigation during the intervention [START_REF] Amiot | Computed tomography-based navigation for hip, knee, andm spine surgery[END_REF]. It offers a unique advantage to visualize the anatomy during surgery and localize anatomical regions without segmenting operative images.

By fusing these images such as CT, C-arm CT, MR or ultrasound with an accurate preoperative model, the surgeon can see the position and orientation of the instrumentation tools on precise anatomical models. While these emerging technologies were successful for knee or hip replacement applications [START_REF] Widmer | Joint replacement-total hip replacement with CT-based navigation[END_REF], corrective spinal surgery is particularly challenging due the complex three-dimensional (3D) deformations of the spine combined with asymmetric deformation of the vertebrae, high variability of the articulated structure and required precision for pedicle screw insertion [START_REF] Kim | Evaluation of pedicle screw placement in the deformed spine using intraoperative plain radiographs with CT[END_REF][START_REF] Lee | Thoracic pedicle screw insertion in scoliosis using posteroanterior C-arm rotation method[END_REF]. Registration of intraoperative fluoroscopic images and preoperative CT/MR images has been proposed to aid interventional and surgical orthopedic procedures [START_REF] Foley | Virtual fluoroscopy: computerassisted fluoroscopic navigation[END_REF]. In some cases, 3D models were registered to 2D X-ray and fluoroscopic images using gradient amplitudes for optimizing the correspondence of single bone structures [START_REF] Markelj | Robust gradientbased 3-D/2-D registration of CT and MR to X-ray images[END_REF]. Objective functions using surface normals from statistical PDMs were applied for the femur [START_REF] Zheng | Unsupervised reconstruction of a patient-specific surface model of a proximal femur from calibrated fluoroscopic images[END_REF] or pelvis [START_REF] Zheng | Statistically Deformable 2D/3D Registration for Accurate Determination of post-operative cup orientation from single standard X-ray radiograph[END_REF]. In spine registration however, one important drawback is that each vertebra of the spine is treated individually instead of as a global shape which hinders surgeons to exploit virtual fluoroscopy imaging. Hence while the morphology of each vertebra remains identical between initial exam and surgery, intervertebral orientation and translation vary substantially.

To tackle this issue, articulated models allow to account for the global geometrical representation [START_REF] Boisvert | Geometric variability of the scoliotic spine using statistics on articulated shape models[END_REF] by incorporating knowledge-based intervertebral constraints. These 3D intervertebral transformations were transposed in [START_REF] Klinder | Spine segmentation using articulated shape models[END_REF] to accomplish the segmentation of the spinal cord from CT images, but multi-modal registration has yet to be solved. Optimization is also based on gradient-descent, prone to non-linearity and local minimums. These methods require segmentation of 3D data or fluoroscopic image, which itself is a challenging problem and has a direct impact on registration accuracy. In [START_REF] Kadoury | Surface/volume-based articulated 3D spine inference through Markov Random Fields[END_REF], we propose a discrete optimization method for articulated spine models, that effectively infers a preoperative model to intraoperative data using only singleton data terms and pairwise constraints. Although intervertebral orientations and translations are nicely captured, the method fails to encode high level geometrical representation of spine changes between the initial exam and during surgery.

In this work, we propose a novel framework which incorporates statistics of regional spinal curves represented by higher order cliques for registering preoperative 3D articulated spine models in a standing position to lying intraoperative 3D CT images. We use a personalized 3D spine reconstructed from biplanar X-rays to derive an articulated model represented with intervertebral transformations. Inference is achieved through a Markov Random Field (MRF) graph which incorporates three optimality components: modular image data-terms which avoids image segmentation, dual potentials to constrain the adjustment of intervertebral links between neighboring objects and global priors of the regional spine components encoded as higher order functionals which are solved INRIA with a primal-dual discrete optimization scheme. We therefore infer a new annotated spine representation from data gathered during surgery based on a high-resolution personalized 3D model.

RT n°0374

2 Personalized 3D Reconstruction of Articulated Spines

Preoperative Spine 3D Reconstruction

From calibrated coronal and sagittal X-ray images I i={1,2} of the patient's spine, the personalized 3D model is achieved by means of a reconstruction method merging statistical and image-based models based on the works of [START_REF] Kadoury | Personalized X-ray 3D reconstruction of the scoliotic spine from statistical and image models[END_REF], and summarized in Fig. 1. The 3D spine centerline C i (u), obtained from quadratic curves extracted from the images is first embedded onto a non-linear manifold containing 732 scoliotic spines (M) to predict an initial spine, modeled by 17 vertebrae (12 thoracic, 5 lumbar), 6 points per vertebra (4 pedicle tips and 2 endplate midpoints). This manifold establishes the patterns of legal variations of spine shape changes in a low-dimensional sub-space based on locally linear embeddings as illustrated in Fig 2 . To map the high-dimensional 3D curve assumed to lie on a non-linear manifold into a low-dimensional subspace, we first determine the manifold reconstruction weights W to reconstruct point i from it's K neighbors, and then determine the global internal coordinates of Y by solving:

Φ(Y ) = M i=1 Y i - K j=1 W ij Y j 2 . (1)
The projection point Y new is used to generate an appropriately scaled model from an analytical method based on nonlinear regression using a Radial Basis Function kernel function f to perform the inverse mapping such that S = [f 1 (Y new ), ..., f D (Y new )] with S = (s 1 , s 2 , . . . , s 17 ), where s i is a vertebra model defined by s i = (p 1 , p 2 , ..., p 6 ), and p i ∈ 3 is a 3D vertebral landmark.

This crude statistical 3D model is refined with an individual scoliotic vertebra segmentation approach by extending 2D geodesic active regions in 3D, in order to evolve prior deformable 3D surfaces by level sets optimization. An atlas of vertebral meshes S i = {x i1 , ..., x iN } with triangles x j are initially positioned and oriented from their respective 6 precise landmarks s i composing S. The surface evolution is then regulated by the gradient map and image intensity distributions [START_REF] Paragios | Geodesic Active Regions: New paradigm to deal with frame partition problems in computer vision[END_REF], where E RAG = αE CAG (S) + (1 -α)E R (S) is the energy function with the edge and region-based components controlled by α are defined as:

E CAG = 2 i=1 Si 1 1 + |∇I i (u i )| α du i (2) E R = - 2 i=1 Πi(Si) log(p R (I i (u i )))du i (3)
with Π i as the perspective projection parameters, and p R is a Gaussian distribution. The projected silhouettes of the morphed 3D models would therefore match the 2D information on the biplanar X-rays in the image domain u, replicating the specifics of a particular scoliotic deformity. At the end of process, the 3D landmark coordinates s i and corresponding polygonal vertebral meshes S i are optimal with regards to statistical distribution and image correspondences.

Articulated Spine Model

The 3D landmarks s i are used to rigidly register each vertebra to its upper neighbor, and the resulting rigid transforms are optimized in the registration problem. Hence, the spine is represented by a vector of local intervertebral rigid transformations A = [T 1 , T 2 , . . . , T N ]. To perform global anatomical modeling of the spine, we convert A into an absolute representation:

A abs = [T 1 , T 1 • T 2 , . . . , T 1 • T 2 • . . . • T N ] (4)
using recursive compositions. The transformations are expressed in the local coordinate system of the lower vertebra, defined by vectors v x , v z and v y = v x × v z , where v x and v z are the vectors linking pedicle and endplate midpoints respectively. Center of transformation is located at the midpoint of all 4 pedicle tips. The rigid transformations described in this paper are the combination of a rotation matrix R and a translation vector t. We formulate the rigid transformation T = {R, t} of a vertebral mesh triangle as y = Rx + t where x, y, t ∈ 3 . Composition is given by

T 1 •T 2 = {R 1 R 2 , R 1 t 2 +t 1 }, while inversion as T -1 = {R T , -R T t}. RT n°0374

Intraoperative Spine Inference from Images with MRFs

A successful inference between the spine model S controlled by the articulations (denoted as A abs ) and the image I must be accomplished by establishing similarity criterions which will drive the model deformation towards the optimal solution. We search the optimal displacement points D = ( d 1 , ..., d n ) of the articulation vectors T that give a good compromise between the encoded prior constraints established by manifold statistics and the fidelity to the image information. Formally, the inference of the model A abs to the image I is given by:

( d 1 , . . . , d n ) = argmin di E(S 0 , I, ( d 1 , . . . , d n )). (5)
The energy E of inferring the spine model S in the image I is a function of the displacement vectors D = ( d 1 , . . . , d n ) in the transformation space applied to the articulation vector A abs . This influences the data-related term V (A 0 abs + D, I) expressing the image cost, a local prior term V (N, D) measuring deformation between neighboring vertebrae and a global higher order term V (H, D) which models the global deformation of a regional curve. The energy function E is therefore:

E S 0 , I, D = V A 0 abs + D, I + V N, D + V H, D = i∈G V i (T 0 i + d i , I) + λ ij i∈G j∈N (i) V ij (T 0 i + d i , T 0 j + d j ) (6)
+ α c c∈C V c (T 0 c + d c )
where A 0 abs + D = {T 0 1 + d 1 , . . . , T 0 n + d n } are the equivalence for articulated components. The data term in the image domain seeks to minimize the distance between model and I:

i∈G V i (T 0 i + d i , I) = Ω η X (I, S i (T 0 i + d i ))dT (7) with η X = vij ∈Si (γ 2 + γ ∇I(v ij ) )/(γ 2 + ∇I(v ij ) 2
) attracts mesh triangles to target high-intensity voxels in the gradient CT volume without segmentation. The term γ is defined as a dampening factor. The second term of Eq.( 6) is a local prior term which are pairwise potentials representing the smoothness term between two consecutive vertebrae and help to constrain the vertebrae main direction in the optimization step (Fig. 3(a)) by assigning a binary clique value V ij = {0, 1} such as:

V ij = |(T 0 i + d i ) -(T 0 j + d j )| ≤ 1 for (i) ≡ (j) || d i || -|| d j || ≤ 2
for (i) = (j). ( 8) The final term in the energy function represents the higher order potentials. We parameterize the potentials with clique variables T c taking on corresponding costs θ q if the cliques are assigned to the displacement vectors d c . To encode anatomical coherence of the spine, three cliques representing each of the spine's regions are composed of three vertebrae, linking the inflexion (change of curvature) and the apical (most deviated) vertebrae as illustrated in Fig. 3(b)-(c). The potential functions are defined as:

INRIA (a) (b) (c)
V c (T 0 c ) = min{ min q∈{1,2,...,t}

θ q + ∆ q (T 0 c ), θ max } (9)
where θ q = ||ψ(T c ) -Π φ(M) (ψ(T c ))|| is a geodesic distance calculated from the low-dimensional mapping of the clique variable onto the manifold using ψ :

→ M to the projection point Π on the prior distribution φ(M), and ∆ q (T 0 c ) = i∈c w q il δ(T i = l) is a deviation function. The Kronecker delta function δ generates binary variables, while the weights are assigned such that:

w q il = 0 if T q = l θ max otherwise ( 10
)
depending on wether the clique variable is given the appropriate label (see sec. 2.3). We develop an optimization procedure which minimizes Eq.( 6) by means of an efficient discrete optimization algorithm using an MRF which is explained in the next section.

Energy Minimization

The optimization strategy for Eq.( 6) of the resulting MRF is based on a discrete labeling principle where we seek to assign the optimal labels L = {l 1 , . . . , l i }, defined in the quantized space Θ = { d 1 , ..., d i } of displacements, to the vertebral transformations represented by nodes T i so that the total energy of the graph is minimum. If we consider that displacing an intervertebral transformation vector by d n is equivalent to assigning label l and that the current solution is given by T t i = T 0 i + t d lit , which adopts the pyramidal coarse-to-fine quantization approach in a temporal minimization, the energy Eq.( 6) can be re-written as a labeling problem:

E t (l 1 , . . . , l n ) = i∈G V i (T t-1 i , l i ) + λ ij i∈G j∈N (i) V ij (T t-1 i , T t-1 j , l i , l j ) (11)
+ α c Tc∈C V c (T t-1 c , l c ).
We solve the minimization of the higher order cliques in Eq.( 11) by transforming them into quadratic functions [START_REF] Rother | Minimizing sparse higher order energy functions of discrete variables[END_REF] using a (t+1)-state switching variable which finds the deviation function which assigns the lowest cost to the labeling:

min V c (T 0 c ) = min T 0 c ,z∈{1,2,...,t+1} f (z) + i∈c g(z, T i ) (12)
where f (z) = {θ q , θ max } is a cost assigning function depending on the state variable z and g(z, T i ) = w q il when z = q and T i = l ∈ L, while g(z, T i ) = 0 when z = t + 1. We apply a Primal-Dual algorithm called FastPD [START_REF] Komodakis | Performance vs computational efficiency for optimizing single and dynamic MRFs: Setting the state of the art with primal-dual strategies[END_REF] which can efficiently solve the registration problem in a discrete domain by formulating the duality theory in linear programming. The advantage of such an approach lies in its generality, efficient computational speed, and guarantees the global optimum without the condition of linearity.

INRIA 5 Experimental Validation

We experimented the articulated inferences of preoperative spine models to intraoperative data by confronting the obtained registered landmark accuracy to expert identification. A dataset of 12 separate CT volumes of the lumbar and main thoracic regions were obtained from different patients (512 × 512 × 251, resolution: 0.8 × 0.8 mm, thickness: 1 -2 mm), acquired for operative planing purposes. Preoperative X-rays of patients were obtained for initial 3D reconstruction. The CT data was manually annotated with 3D landmarks, corresponding to left and right pedicle tips as well as midpoints of the vertebral body. The smoothness term was set at λ ij = 0.4, while the clique variable was α c = 0.3. Tests were performed in C++ on a 2.8 GHz Intel P4 processor and 4 GB memory.

We first experimented the manifold-based higher order energy functions based on the search of the label space. Since the functional potentials are directly linked to the costs assigned by the geodesic distances to the manifold, we can therefore assess the performance of this metric as an efficient parametrization of the higher order cliques. Fig. 4 displays the obtained embedding for a dataset of 711 spine models in M, as well as selected energy potential distributions. Results show the search for the minimum energy is a able to obtain cliques that fall near the manifold, thus corresponding anatomically coherent configurations.

Then for each case of the CT dataset, registration is performed to automatically align the CT volume with γ = 0.05 to the given preoperative model S and quantitative assessment consisted of measuring the RMS distance with the manually segmented landmarks. Table 1 presents the results from this experiment with 3D landmark RMS differences and surface errors based on DICE scores, comparing the proposed method with cost functions excluding pairwise and higher order energy terms. Results for overall vertebral landmark errors have improved by 0.38 mm compared to the previous approach [START_REF] Kadoury | Surface/volume-based articulated 3D spine inference through Markov Random Fields[END_REF], which is significant for the required guidance accuracy. These results seem to confirm that exploiting higher order geometrical constraints on the whole shape prior does help to converge towards a global minimum. Visual registration results of the 3D model with CT is shown in Fig. 5, demonstrating the multi-modal alignment where one could observe accurate superposition of geometrical models on selected multi-planar views. Table 1: Comparison between optimization schemes with cost functions integrating only image related terms, singleton-pairwise potentials, and the proposed method with additional higher order functions. Evaluation is made on root-mean-square (RMS) differences (inferred vs. annotated landmarks) and with DICE scores from segmented vertebrae.

Conclusion

Statistical deformable model methods are often dedicated to single anatomical structures. Shape analysis of articulated models on the other hand has been sparsely investigated due to the difficulty in constraining the higher number of transformation variables. The method we propose not only allows to infer shape deformations of object constellations using discrete optimization techniques, but offers the possibility to learn the variations of spinal shape in complex corrective procedures. Hence we progress to whole body deformation by introducing novel high order cliques in the optimization to impose global shape constraints.

Increased accuracy of pedicle landmark localization was achieved by comparing the method with standard image-based methods. Future work will look at adapting local shape variations using mesh relaxation techniques based on bone density fields which may increase the accuracy of the geometrical alignment. Furthermore, a real-time feasibility evaluation of the approach during corrective surgery is planned for clinical use.

RT n°0374

INRIA

Figure 1 :

 1 Figure 1: Personalized spine 3D reconstruction from preoperative X-rays.

Figure 2 :

 2 Figure 2: Illustration of spine distribution embedded onto on a low-dimensional manifold.

Figure 3 :

 3 Figure 3: (a) Pairwise link between neighboring vertebrae; (b) Regional higher order cliques based on (c) vertebral triplets.

Figure 4 :

 4 Figure 4: (a) Low-dimensional manifold of the spine dataset used to determine costs for clique variables. (b) High order energy potential distributions. (c) Parametrization of manifold clique variables in the label search space corresponding to the original energy function.

Figure 5 :

 5 Figure 5: Qualitative assessment of multimodal registration results of spine model S with CT images I.

  Inférence de Modèles 3D Préopératoire de la Colonne Vertébrale aux Données Volumétriques par des MRFs de Haut-Niveau Résumé : Ce papier présente une méthode d'inférence d'un modèle personnalisé de la colonne vertébrale en 3D à partir de données tomodensitométriques en exploitant des fonctions d'énergie de haut niveau incorporant des dépendances géométriques. Une reconstruction 3D précise à partir d'images radiographiques standards est exploitée afin d'obtenir une représentation du rachis modélisée par une série de transformations intervertébrale. Ces transformations sont basées sur les paramtres de rotation et de translation. La transformation de la forme du rachis entre les positions couchées et debout est atteinte grâce a une optimisation un Markov Random Field (MRF), où les variables inconnues sont les déformations appliquées aux transformations intervertébrales. Des valeurs potentiels unitaires et binômes mesurent le lien entre les images et les contraintes géométriques entre les vertèbres, alors que des fonctions de haut niveau introduisent des contraintes d'alignement des régions anatomiques.

L'optimisation des paramètres dans un contexte multi-modale est effectuée par une approche de programmation linéaire et par dualité. Nous présentons des résultats prometteurs pour le recalage d'images à partir d'une comparaison avec une identification manuelle d'un expert qui offre une précision de 1.8 ± 0.7mm basée sur la localisation de repères chirurgicaux.

RT n°0374

http://www.inria.fr ISSN 0249-0803


Martin Deraux
email: deraux@ujf-grenoble.fr

Harish Seshadri
email: harish@math.iisc.ernet.in

ALMOST QUARTER-PINCHED K ÄHLER METRICS AND CHERN NUMBERS

Given n ∈ Z + and ε > 0, we prove that there exists δ = δ(ε, n) > 0 such that the following holds:

and c I (M ), c J (M ) are any two Chern numbers of M , then

where c 0 I , c 0 J are the corresponding characteristic numbers of a complexhyperbolic space-form.

It follows that the Mostow-Siu surfaces and the threefolds of Deraux do not admit Kähler metrics with pinching close to 1 4 .

Introduction

It is well known that Riemannian manifolds of constant sectional curvature -1 must be locally isometric to real hyperbolic space. Allowing the sectional curvatures to be merely close to -1 (say with pinching close to 1, i.e. with a suitable normalization all sectional curvatures are between -α and -1 for some α close to 1), there are many more examples than hyperbolic manifolds.

In fact, in [START_REF] Gromov | Pinching constants for hyperbolic manifolds[END_REF], Gromov and Thurston construct manifolds of dimension ≥ 4 that admit negatively curved Riemannian metrics with pinching arbitrarily close to 1, but that do not admit any metric of constant curvature (in fact they do not admit any locally symmetric metric).

The Kähler analogue of the above examples is not well understood. Negatively curved Kähler metrics must have pinching at most 1 4 , as can be deduced from an inequality due to Berger (see [START_REF] Berger | Sur quelques variétés Riemanniennes suffisamment pincées[END_REF]), and 1 4 -pinched metrics are precisely those locally isometric to complex hyperbolic space (of dimension ≥ 2). It is not known, however, whether there exist negatively curved Kähler manifolds with pinching arbitrarily close to 1 4 which are not biholomorphic to ball quotients. The surfaces constructed by Mostow and Siu in [START_REF] Mostow | A compact Kähler surface of negative curvature not covered by the ball[END_REF] (as well as their threedimensional analogue, see [START_REF] Deraux | A negatively curved Kähler threefold not covered by the ball[END_REF]) admit Kähler metrics with negative sectional curvatures, but very little is known about the pinching of these metrics.

The main result of this note relates proximity to 1 4 -pinching with proximity of the ratios of Chern numbers to the corresponding ratios of Chern numbers of a complex hyperbolic manifold of the same dimension. This shows that neither the Mostow-Siu surfaces, nor the three-dimensional analogues constructed in [START_REF] Deraux | A negatively curved Kähler threefold not covered by the ball[END_REF], admit Kähler metric with pinching close to 1 4 . This answers a question raised by Pansu

Mathematics Subject Classification (1991): Primary 53C21, Secondary 53C20.

in [START_REF] Pansu | Pincement des variétés à courbure négative d'après M. Gromov[END_REF] (echoing a question of Gromov), which seems not to be answered anywhere in the literature.

Pinching and ratios of Chern numbers

The main observation of this note is the following:

Theorem 2.1. Given n ∈ Z + and ε > 0, there exists δ = δ(ε, n) > 0 such that the following holds:

If (M n , g) is a compact Kähler n-manifold whose sectional curvatures K satisfy -1 -δ ≤ K ≤ - 1 4 ,
and c I (M ), c J (M ) are any two Chern numbers of M , then

c I (M ) c J (M ) - c 0 I c 0 J < ε,
where c 0 I , c 0 J are the corresponding characteristic numbers of a complex-hyperbolic space-form.

Recall that Mostow and Siu constructed an infinite family of surfaces that admit negatively curved Kähler metrics, but the corresponding ratios c 2 1 /c 2 do not approach 3 as one ranges through these examples (see the introduction of [START_REF] Mostow | A compact Kähler surface of negative curvature not covered by the ball[END_REF]). The ratios c 3 1 /c 3 for the threefolds constructed in [START_REF] Deraux | A negatively curved Kähler threefold not covered by the ball[END_REF] certainly cannot approach 16, since only a finite number of threefolds were constructed there.

Theorem 2.1 implies that the Mostow-Siu surfaces (as well as the threefolds in [START_REF] Deraux | A negatively curved Kähler threefold not covered by the ball[END_REF]) do not carry Kähler metrics with pinching close to 1/4.

It is an interesting open question whether the conclusion of Theorem 2.1 holds if only assumes the existence of almost-quarter pinched Riemannian metrics on a compact Kähler manifold. This question is natural in light of the result of Hernandez-Zheng-Yau ([4], [START_REF] Yau | Negatively 1 4 -pinched Riemannian metric on a compact Kähler manifold[END_REF]): If a compact Kähler manifold M admits a negatively quarter-pinched Riemannian g metric, then (M, g) is holomorphically isometric to a complex-hyperbolic space-form.

proof

Let V be a vector space over R of real dimension 2n with an almost complex structure J : V → V and a J invariant inner product , . Hence J 2 = -I and J(v), J(w) = v, w for all v, w ∈ V .

A curvature tensor on V is a 4-linear map R :

V × V × V × V → R, i.e. an element of ⊗ 4 V * , satisfying (1) R(X, Y, Z, W ) = -R(Y, X, Z, W ) = -R(X, Y, W, Z) (2) R(Z, W, X, Y ) = R(X, Y, Z, W ) (3) R(X, Y, Z, W ) + R(X, W, Y, Z) + R(X, Z, W, Y ) = 0. If, in addition, R satisfies (4) R(J(X), J(Y ), Z, W ) = R(X, Y, J(Z), J(W )) = R(X, Y, Z, W )
then we say that R is a Kähler curvature tensor.

Let R ⊂ ⊗ 4 V * denote the space of Kähler curvature tensors and R 0 ∈ R denote the curvature tensor of complex hyperbolic space CH n .

It is clear from this expression and the earlier one expressing K(x, y) in terms of holomorphic sectional curvatures that |R(x, y, z, t) -R 0 (x, y, z, t)| < ε (2n) 4 if |H(u) + 1| < η for an η = η(ε, n). Proposition 3.2. Given ε > 0, there exists δ = δ(ε, n) such that if the sectional curvatures of a Kähler curvature tensor R ∈ R satisfy

-1 -δ ≤ K ≤ - 1 4 then |R -R 0 | < ε in the norm on ⊗ 4 V * .
Proof. By the lemma above, it is enough to show that there is a δ = δ(η) > 0 such that -1δ ≤ K ≤ - 1 4 implies |H(u) + 1| < η for all unit vectors u. Let u be any unit vector in V . Choose v so that {u, J(u), v, J(v)} is an orthonormal set. By ( 3) and ( 4) in the definition of a Kähler curvature tensor,

(3.2) K(u, v) + K(u, J(v)) -R(u, J(u), v, J(v)) = 0.
We recall the following estimate of Berger ([1], Page 69, ( 7)). Note that though Berger works with positive pinching, his proof works for negative pinching as well. Also, our pinching hypothesis and the conclusion can be obtained from his by a rescaling. If R is a curvature tensor (not necessarily Kähler) which is α-pinched,

i.e. -α ≤ K ≤ -1 4 then |R(X, Y, Z, W )| ≤ 2 3 (α - 1 4
).

for any orthonormal set {X, Y, Z, W }.

Let η > 0 be as in Lemma 3.2. By Berger's estimate there exists δ 1 = δ 1 (η) such that if the sectional curvatures K satisfy

(3.3) -1 -δ 1 ≤ K ≤ - 1 4 then (3.4) |R(u, J(u), v, J(v))| ≤ 1 2 + η 6 .
By (3.2) and (3.4)

K(u, v) = -K(u, J(v)) + R(u, J(u), v, J(v)) ≥ 1 4 - 1 2 - η 6 = - 1 4 - η 6 (3.5)
Of course, we have the same estimates for K(u, J(v)), K(J(u), v)) and K(J(u), J(v)). Hence,

(3.6) - 1 4 - η 6 ≤ c ≤ - 1 4 , - 1 2 - η 3 < R(u, J(u), v, J(v)) ≤ - 1 2 ,
where c ∈ {K(u, v), K(u, J(v)), K(J(u), v), K(J(u), J(v))}.

Applying the second inequality in (3.6) to the orthonormal set

{ 1 √ 2 (u -J(v)), 1 √ 2 (J(u) + v), 1 √ 2 (u + J(v)), 1 √ 2 (J(u) -v)}, we get (3.7) R(u -J(v), J(u) + v, u + J(v), J(u) -v) ≤ -2.
Expanding the left hand side

(3.8) R(u -J(v), J(u) + v, u + J(v), J(u) -v) = H(u) + H(v) + 2R(u, J(u), v, J(v)) -4K(u, v) Combining (3.6
) and (3.7) we have

H(u) + H(v) -4K(u, v) < -1 + 2η 3 . ( 3
.9) (3.9) and the first inequality in (3.6) imply

H(u) + H(v) < -2 + 2η 3 .
Finally, using (3.3), the above inequality implies

-1 -δ 1 < H(u) < -1 + 2η 3 + δ 1 . Let δ = min{ η 3 , δ 1 }. Then |H(u) + 1| < η
and the proof is complete since the unit vector u was arbitrary.

Let ω denote the Kähler form of V , defined by ω(u, v) = u, J(v) Then the volume form on V , with V carrying the natural orientation coming from J, is just ω n .

Let I = (a 1 , ..., a n ) and J = (b 1 , ..., b n ) be muti-indices in (Z ≥0 ) n with Σ i ia i = Σ i ib i = n. Given an Kähler curvature tensor R we can form the corresponding Chern forms c I (R) = c a1 1 (R) ∧ ... ∧ c an n (R) and c J (R) = c b1 1 (R) ∧ ... ∧ c bn n (R). These are elements of ∧ n V * .

Given ε 1 (to be chosen below), Proposition (3.2) (and the obvious continuous dependence of the Chern forms on the full curvature tensor) implies there is a δ = δ(ε 1 , n) such that -1δ < K < -1 implies that (3.10)

|c I (R) -c I (R 0 )| < ε 1 , |c J (R) -c J (R 0 )| < ε 1 in the norm on ∧ n V * . Let c I (R 0 ) = a(n)ω n and c J (R 0 ) = b(n)ω n . Choose ε 1 = ε 1 (ε, n) so that (3.11) a(n) -ε 1 b(n) + ε 1 > a(n) b(n) -ε a(n) + ε 1 b(n) -

  n)ε 1 )ω n < c I (R) < (a(n) + ε 1 )ω n and (b(n)ε 1 )ω n < c J (R) < (b(n) + ε 1 )ω n

We now turn to a compact Kähler manifold with almost-quarter pinching as in Theorem 1.1. We can integrate the above inequalities on M to get

V ol(M ). Dividing these two inequalities and using (3.11) gives the desired result.

R 0 is given by -4R 0 (u, v, z, w) = u, z v, wu, w v, z + u, J(z) v, J(w) (3.1) u, J(w) v, J(z) + 2 u, J(v) z, J(w)

It is well-known that if R is a Kähler curvature tensor on V whose sectional curvatures are negative and 1 4 -pinched, then R = R 0 . We generalize this in Proposition 3.2. Recall that we regard the space R of Kähler curvature tensors as a subspace of ⊗ 4 V * . The inner product on V canonically induces an inner product on ⊗ 4 V * and hence on R.

In what follows, we will use the following notation: For arbitrary u, v ∈ V ,

If u and v are orthornormal, this is the sectional curvature of the 2-plane spanned by u and v. For any u, let

If u is a unit vector this is the holomorphic sectional corresponding to u. Note that

for every orthonormal set {x, y, z, t}.

Let u, v ∈ V be unit vectors. If a, b are real numbers such that a 2 + b 2 = 1, then (this calculation is from [START_REF] Kobayashi | Holomorphic bisectional curvature[END_REF], Lemma 3)

We also have the equation

which follows from (3) and (4) in the definition of a Kähler curvature tensor. Take

It is clear that we can use these three equations to express R(u, v, u, v) as a linear combination of

Moreover note that the coefficients in this expression are absolute constants, not dependent on anything. On the other hand, one knows that the full curvature tensor can be expressed in terms of sectional curvatures. In fact [START_REF] Mendonca | Expression of curvature tensors and some applications[END_REF], 24R(x, y,z, t) = K(x + z, y + t) + K(xz, yt) -K(x + z, yt) -K(xz, y + t) -K(x + t, y + z) -K(xt, yz) + K(x + t, yz) + K(xt, y + z)


Abdulhalim Dandoush
email: abdulhalim.dandoush@sophia.inria.fr.

Alain Jean-Marie

Thème

Download Process in Distributed Systems, Flow-level Algorithm vs. Packet-level Simulation Model

Keywords: distributed systems, performance evaluation, download time, response service time, data distribution, simulation model, max-min fairness systèmes distribués, évaluation de performance, temps de téléchargement, distribution de données, modèle de simulation, équité max-min

Parallelism in the download process of large files is an efficient mechanism for distributed systems. In such systems, some peers (clients) exploit the power of parallelism to download blocks of data stored in a distributed way over some other peers (servers). Determining response times in parallel downloading with capacity constraints on both the client downloads and server uploads necessitates understanding the instantaneous shares of the bandwidths of each client/server is devoted to each data transfer flow.

In this report, we explore the practical relevance of the hypothesis that flows share the network bandwidth according to the max-min fairness paradigm. We have implemented into a flow-level simulator a version of the algorithm which calculates such a bandwidth allocation, which we have called the "progressive-filling flow-level algorithm". We have programmed a similar model over NS2 and compared the empirical distributions resulting from both simulations.

Our results indicate that flow-level predictions are very accurate in symmetric networks and good in asymmetric networks. Therefore, PFFLA would be extremely useful to build flow-level simulators and, possibly, to perform probabilistic performance calculations in general P2P networks.

Introduction and related work

The growth of storage volume, bandwidth, and computational resources for PCs has fundamentally changed the way applications are constructed. Almost 10 years ago, a new network paradigm has been proposed where computers or peers can build a virtual network (called overlay) on top of another network or an existing architecture (e.g. Internet). This new network paradigm has been labeled peer-to-peer (P2P) distributed network. A peer in this paradigm is a computer that play the role of both supplier and consumer of resources, in contrast to the traditional client-server model where only servers supply, and computers consume. Applications that use this distributed network provides enhanced scalability and service robustness as all the connected computers or peers provide some services.

This distributed network model has proved to be an alternative to the Client/Server model and a cheap, scalable, self-repairing and promising paradigm for grid computing, grid delivery network (GDN), file sharing, voice over IP (VoIP), backup and storage applications.

Such distributed systems rely on data fragmentation and distributed storage. Files are partitioned into fixed-size blocks that are themselves partitioned into fragments. Fragments are usually stored on different peers. Given this configuration, a user wishing to retrieve a given block of data would need to perform multiple downloads, generally in parallel for an enhanced service. The transfer of sequences of packets on one long-term TCP connection (e.g. download a fragment of data between two peers in a P2P system or between a client and a server through FTP protocol) defines a "flow". A flow can as well refer to the sequences of packets that constitute a block of data and that follow several TCP connections simultaneously. In this work, we will consider the former definition.

One measure of the quality of the service given by the distributed storage/parallel download infrastructure is the time it takes to retreive the complete document. This in turns depends on the throughput of the different flows created to obtain the fragments of this document. Their values are, a priori, a function of the demand and capacities of the complete network entities: clients, servers and links.

The basic problem of predicting the instantaneous shares of the bandwidth received by each flow of a TCP-based network has received quite some attention in the last 15 years, in connection with the notion of fairness; yet, there is no clear consensus in the literature on a simple formula or algorithm to give a reasonable solution of this problem. Such an algorithm would be extremely useful to build flow-level simulators and, possibly, to perform probabilistic performance calculations.

On the one hand, some authors have shown that the dynamics of TCP have been shown to be quite chaotic is some situations. Other authors on the other hand, have argued that TCP tends to share the bandwidth between flows reasonably. For instance, Heyman et al. [START_REF] Heyman | A new method for analysing feedbackbased protocols with applications to engineering web traffic over the internet[END_REF], followed by Fredj et al. [START_REF] Ben Fredj | Statistical bandwidth sharing: a study of congestion at flow level[END_REF], have studied a single bottleneck link shared by a given number of identical sources that alternately send documents though the shared link and stop sending for randomly thinking time. They showed through simulations that TCP shares fairly the bottleneck (that is, in equal shares) and they introduced analytical tools that can predict the expectation of the transmitting rate. Varki proposed in [START_REF] Varki | Response time analysis of parallel computer and storage systems[END_REF] a simple approximation for the expected response time based on the fork-join model. Massoulié and Roberts proposed in [START_REF] Massoulié | Bandwidth sharing and admission control for elastic traffic[END_REF] a model similar to that of [START_REF] Heyman | A new method for analysing feedbackbased protocols with applications to engineering web traffic over the internet[END_REF] where the interflows arrival times are i.i.d. exponentially distributed random variables. They studied the network as M/G/1 PS queue. In [START_REF] Chiu | Minimizing file download time in stochastic peer-to-peer networks[END_REF], Chiu and Eun, the authors have focused on the average download time of each user in a P2P network while considering the heterogeneity of service capacities of peers. They point out that the common approach of analyzing the average download time based on average service capacity is fundamentally flawed.

Other studies have put forward the concepts of max-min fairness, proportional fairness, balanced fairness and utility-based resource-sharing models (see e.g. [START_REF] Boudec | Rate adaptation, Congestion Control and Fairness: A Tutorial[END_REF] and the reference therein). One conclusion of these studies is that throughput allocations resulting from the use of the TCP protocol for infinitely long flows are usually not max-min fair. However, the results of Bonald and Proutière [START_REF] Bonald | Insensitive bandwidth sharing in data networks[END_REF] suggested that when the flows are dynamic (flows are continuously created and have a finite duration), the average throughput obtained by flows under various sharing mechanism tend to be similar. It is quite possible that, from a practical perspective, the predictions obtained with a maxmin fair sharing mechanism may be "good enough".

The purpose of this paper is to assess whether max-min fairness for the allocation throughput is a proper model when evaluating response times of parallel downloads.

Given the variety of situations to be studied, we begin with the simplest scenario: a symmetric network in which we assume that capacity constraints are located at the client/server nodes, and not inside the network. We also assume that all RTTs are equal. 1We use an algorithm which calculates an instantaneous throughput for each individual flows in a certain set of flows, given the upload and download capacities of the client and server nodes. This algorithm can be seen as a variant of the "progressive filling" [START_REF] Boudec | Rate adaptation, Congestion Control and Fairness: A Tutorial[END_REF] or "water filling" algorithm of [START_REF] Bertsekas | Data Networks[END_REF]: we name it as the Progressive Filling Flow Level Algorithm (PFFLA). The validation of this algorithm consists in characterizing the response time of parallel downloads in a distributed storage system, through simulations. We have implemented the PFFLA in a flow-level simulator of parallel downloads, and we have programmed a similar model over NS2. The response times in the flow level simulator have been compared to that of the packet-level simulations in NS (both distributions and averages). This experimental setting is, to the best of our knowledge, original in at least three features. First, we consider flows related to downloads in parallel, which are synchronized when they are created. Second, we consider that the possible bottlenecks for flows occur only at the edge of the network, never inside. Finally, we consider large numbers of nodes (up to 500) and flows.

Our results show that the relative error between PFFLA and NS-2 for the expected value is less than 2% for relatively large loads in the system (e.g. ρ = 60%) and less than 1% for low loads in the symmetric case and less than 7% respectively for relatively large loads in the system (e.g. ρ = 50%) and less than 1% for low loads in the asymmetric case. We conclude that PFFLA is a reliable mechanism to analyze the service response time in many systems based on P2P and Grid computing concepts such as Storage Systems and Grid Delivery Networks.

The rest of this paper is organized as follows. Section 2 overviews the system assumptions and notation. Section 3 describes the flow-level simulation algorithm "PFFLA". In Section 4, comparisons between packet-level and PFFLA are introduced and discussed. Last, Section 5 concludes the paper and highlights some future directions.

System description and notation

In the following, we will distinguish the servers, which are computers that provide a storage service, from the clients whose objective is to retrieve data from the servers to account for the fact that flows (transfer of sequences of data units from a server to a client) have a direction. In the terminology of P2P-based systems, each "peer" has the role of both a client and a server. It is usual that the communication link from the network to the peer (upload link) and the one from the peer to the network (download link) are not shared. Their capacity may actually not even be the same, as with ADSL network accesses. In that case, the entities client and server can be considered as two distinct nodes. On the other hand, if the network access is indeed shared between input and output, the peer is represented by one node. In the following, we shall only consider the first situation.

In this study, we are interested in systems where blocks of data are partitioned into several equally sized fragments stored randomly over different servers. We will consider both homogeneous (upload/download capacities are identical) and asymmetric situations.

We consider a distributed system in which the following assumptions and notations will be enforced throughout the paper.

Network assumptions

• The considered network consists in a set of nodes N . In a P2P context, there will be N /2 peers according to this notation. In other words, we will have N /2 servers and N /2 clients.

• The logical structure of the network is that of a star, with an infinite-capacity central node. In other word, the interconnection network underlying the parallel download application is assumed not to introduce capacity constraints. Only the upload or download links (the branches of the star) have a limited capacity.

The capacity of upload links (from servers to the network) is C u , the capacity of download links (from network to clients) is C d .

• The temporal distance (measured as the round-trip time, RTT) is assumed to be the same between pairs of nodes (clients or server).

Data and traffic assumptions

• Each block of data D of size S B is partitioned into s fragments of size S F .

• We assume that the s servers that hold fragments of a given block of data D are uniformly selected over all servers in the system, and are all distinct.

• Each download request of a block of data issued by a client will generate s parallel requests toward s servers to retrieve s distinct fragments of the requested block. A request generates s flows.

The assumption on the uniform distribution of the blocks of some document corresponds to the situation where a very large number of documents exist, and/or each fragment of each document has been replicated a large number of times. In that situation, it is unlikely that the set of blocks needed by two distinct requests will be correlated. The network being symmetric, it is reasonable to assume that fragments have been uniformy distributed. The assumption that different fragments of some document are stored on different peers is common in P2P-based systems: it results mainly from privacy and data ownership issues.

Description of the algorithm

Before describing the algorithm we have used, we recall the principle of max-min fairness, and the algorithm (referred to as the "Progressive Filling" algorithm in [START_REF] Boudec | Rate adaptation, Congestion Control and Fairness: A Tutorial[END_REF]).

The notion of max-min (or maximin) fairness originates from the field of political philosophy and economics, and was introduced in the context of networking by Bertsekas and Gallager [2, ch. 6] as a design objective for communication networks, in particular, the design of flow control schemes. The main idea of max-min fairness is to maximize the allocation of each flow f subject to the constraint that an incremental increase in f 's allocation does not cause a decrease in some other flow's allocation that is already as small as f 's or smaller [2, p. 526].

For the purpose of formalizing the description of algorithms, introduce the following notation. The network is assumed to be made of a set A of links. Each link a has a capacity C a . The traffic is formed by a set F of flows. We assume that flows cannot be split between several routes of the network. This implies that we can assume that each flow f has a throughput θ f ≥ 0, and crosses certain links of A. We write f ∇a to denote the fact that f crosses a.

Using this notation, the total flow on link a of the network is then given by:

F a = f ∇a θ f .
The capacity constraint for the network is then:

F a ≤ C a , ∀a ∈ A . (1)
A vector of throughputs {θ f } satisfying these constraints is said to be feasible.

The existence of a max-min fair allocation of throughputs is not entirely obvious. Indeed, satisfying the capacity constraints implies that the increase of some flow's throughput may result in the decrease of another flow's throughput, and conversely.

Bertsekas and Gallager have proved that there exists one unique feasible allocation of throughputs which is max-min fair. It can be constructed using the following algorithm.

The idea is to start with an all-zero rate vector and to increase the rates on all paths together until F a = C a for one or more links a. At this point, each flow using a saturated link has the same throughput at every other flow using that link. Thus, these saturated links serve as bottleneck links for all flows using them.

INRIA

At the next step, all flows not using the saturated links are incremented equally in rate until one or more new links become saturated. The newly saturated links serve as bottleneck links for those flows that pass through them but do not use the previously saturated links. The algorithm continues until all flows pass through at least one saturated link. This process is often visualized as progressively augmenting the throughput until capacities are "filled", hence the name of "progressive filling". Information flows are also sometimes visualized as some "fluid" which is poured into the network. For this reason, the algorithm is also referred to as a "water filling" algorithm.

We have used in our analysis the following version of this algorithm.

Data: A set of links A with their capacities C a , and a set of flows F Result: A throughput value for each flow, satisfying the throughput constraint (1) begin

Remove from A nodes without flows ;

while A not empty do foreach a ∈ A do N a ← #{f ∈ F |f ∇a} ; end calculate θ * = min a∈A C a /N a ; calculate a * = arg min a∈A C a /N a ; foreach f , f ∇a * do set θ f = θ * ; foreach a, f ∇a do C a ← C a -θ * end remove f from F ; end remove from A links without flows ; end return {θ f } end Algorithm 1: Algorithm PFFLA
The fact that this algorithm produces a max-min allocation can be checked the same way as for the Progressive Filling Algorithm: according to [2, p. 527], max-min solutions are characterized by the "bottleneck" property: ∀f ∈ F, ∃a ∈ A, f ∇a and g∈F ,g∇a

θ g = C a and ∀h ∈ F, h∇a, θ h ≥ θ f . (2)
This property can be checked almost by construction on our algorithm.

More remarks concerning this algorithm:

• The algorithm eventually stops because at least one link is removed from A at each loop. The number of loops is therefore bounded by the cardinal of the initial set of links. Since several links can be removed in each loop, the algorithm may actually stop faster.

• When arg min a∈A C a /N a contains more than one element, it does not matter which one is chosen, in the sense that the outcome of the algorithm does not depend on that particular choice.

This results, by contradiction, from the fact that the max-min allocation is unique. It can also be proved that other links of the set are still in the "argmin" at the next step, so that each of them will be chosen eventually. Equivalently, one may remove simultaneously from the network all flows that are connected to links in the "argmin" set.

• It is possible to add constraints on the throughput of flows. For instance, the throughput of a TCP flow on a lossless connexion with RTT τ and maximum window size w is always less than w/τ .

In our situation, we have made the assumption that the network can be represented by a star, and the flows cross exactly two links: one upload link and one download link. The algorithm is capable to handle general situations however.

Experimental results

Parameter values

We ran a total of fifteen experiments; ten in symmetric peers download/upload capacities scenarios (homogeneous networks) and five in asymmetric scenarios.

The set-up of the common parameters between the two simulation levels is summarized in Table 1. The capacities that we have selected in the simulations vary between the values of the ISDN and ADSL technologies (384, 576 and 1500 kbps). In experiments 1-10, nodes are homogeneous: they have all the same network access capacity. In Experiments 11-17, capacities of clients and servers are asymmetric. Download requests at each client node arrive according to some Poisson process of given rate λ. The different request processes are independent. This assumption is reasonable in practice: Guha et al. have shown in [START_REF] Guha | An experimental study of the skype peer-to-peer VoIP system[END_REF] that in real networks, and when the number of clients N c is large, the request arrival process can be reasonably modeled by a Poisson process. We vary the value of the request generation rate across the experiments such that the total load in the system ρ (see below) varies from low (e.g. 6%) up to high value (e.g. 70%) as reported in Table 1.

The last setting concerns the blocks and fragments sizes that are stored in the system. Fragment sizes S F (resp. block sizes S B ) in P2P systems, for instance, are typically between 256KB and 4MB each (resp. between 4MB and 9MB each). We will consider in most of our experiments S F = 2MB and S B = 8MB, except in Experiment 1 where S F = 1MB and S B = 4MB. Therefore s = 4 in all experiments. In the asymmetric scenarios, we have chosen the two values 1500/384 so that the capacity of a server is slightly larger than 1/s times the capacity of a client.

INRIA

For the packet-level simulation details, we consider a fixed constant value of 2ms for the link propagation delays. The main TCP configurations are as follow: we use TCP segment size (S pkt ) of 1460 Bytes, the upper bound on the advertised window for the TCP connection is set to 40, the initial size of the congestion window on slow-start is 2, and the TCP/IP header size (h ip ) is 40 Bytes. The P2P application layer header (h a ), which is implemented over the NS transport layer, is 13 Bytes for each fragment. The queue management type used in the links is "DropTail" with size of 500 packets. The maximum window size is left to NS2's default of 64kB. Given our assumptions on propagation delay, this gives a maximum TCP throughput of 64kB/8ms = 4M B/s, largely superior to the capacity of the links. Therefore, maximum window effects are not expected to restrict the throughput of file transfers.

In the flow-level simulation, and when calculating the total amount of data sent in the TCP flows, we neglect the fact that one data packet may be incomplete after segmentation. We also neglect the packets sent during the opening and the closing of the TCP connection, and we assume that no retransmission occurs. The total amount of data transported during the download of one document is then calculated by multiplying the application data size by the overhead factor due to packet headers, that is:

L(bits) = s × (S F (bits) + h a (bits)) × (1 + h ip /S pkt ) . (3)
Consider a client node with link capacity C. The time to download a complete document would be, when no interferences from other downloads occur:

σ = s × (S F (bits) + h a (bits)) × (1 + h ip /S pkt ) C(bps) . (4)
On the other hand, if the global arrival rate of document requests is λ, the rate of requests arriving at a particular client is λ/(N /2). Accordingly, the load factor of a client link of capacity C in the network is:

ρ = λs × (S F (bits) + h a (bits)) × (1 + h ip /S pkt ) C(bps)N /2 . ( 5
)
Consider now a server node with link capacity C. Given our assumption on the uniform repartition of blocks on servers, the rate of arrivals of fragment requests at the servers is λs/(N /2). The duration of one request should be σ/s since only one fragment is concerned. Finally, the load factor of the server's link is given by Equation ( 5) also.

In the homogeneous cases, this value of ρ can also be interpreted as the load of the whole network (ratio of global data requests to global transfer capacity). In the asymmetric cases, we take ρ as the load of the links with the smallest capacity.

Simulators and Metrics

We have developed a packet-level simulator and a flow-level simulator for our model. The packetlevel simulator is build using NS2. Its implementation details can be found in [START_REF] Dandoush | A realistic simulation model for peer-to-peer storage systems[END_REF].

The flow-level simulator consists in the embedding of Algorithm 1 into a discrete-event simulator handling the arrival and the departure of flows. The principle is that every time the set of flows present in the network changes, the bandwidth shares are re-computed with the algorithm, and it is assumed that these throughputs are obtained instantaneously. The program keeps track of the remaining quantities to be downloaded in each flow, and can compute the date of the next event: arrival or end of download.

Both simulators are instrumented so as to produce response times for fragments and complete documents.

The metric we are interested in is the download time of a document. For a given request, this is the maximum between the download time of the s fragments of the document. Of course, this is a random variable, and we measure its empirical distribution and empirical average.

In addition to the simulations, we have compared the average document download time, denoted by E[T N S ] which is calculated from the packet-level simulation, with the average response time in a simple queueing system. The rationale for this is that, if the throughput of the connections were limited only by the client's capacity, then the link would behave as a Processor Sharing queue. This is because the size of the fragments is the same, so that the response time of all s fragments is the same, and all s fragments can be actually considered as a single "client". The client's bandwidth is then shared between different requests. Since requests arrive according to a Poisson process, the model is that of a M/D/1 processor sharing (PS) queue. This model is expected to work well when the load is small: indeed, in that case it is unlikely that flows will be limited on the server side. Since the average response time in this queue, denoted by E[T P S ], can be computed with a simple formula, we can easily test this conjecture.

The average response time in the M/D/1/P S queue is known to be (see e.g. [START_REF] Kleinrock | Queueing Systems[END_REF]):

E[T P S ] = σ 1 -ρ = s(•S F (bits) + h a (bits)) × (1 + h ip /S pkt ) C(bps) × (1 -ρ) (s) . (6)
We will compute the relative error (RR) between E[T N S ] and E[T P S ], on one hand, and between E[T N S ] and E[T F LA ] on the other hand.

Results on the PS model also include the distribution of the response time. The relevant formulas are provided in the Appendix.

Results

We have run both the flow-level simulator and the packet-level simulator on the seventeen sets of values described in Table 1.

For the flow-level simulations, we have collected 100000 samples of the document download time in every case. The number of samples collected with the packet-level simulation is reported in Table 1 varies from 30000 to 70000 samples for most experiments. In Experiment 10, we collected roughly 280000 samples.

We also report in Table 1 the execution time of the packet-level simulations. The experiments were run on a machine with the following principal characteristics: multithreaded processor Intel(R) Core(TM)2 Duo of 2.66GHz, 4GB RAM + 4GB swap running Fedora Core 5. The analysis of number of samples issued by unit time of computation (figure not reported) reveals that this number decreases with the number of nodes. Interestingly, it tends to slowly decrease when the load ρ of the INRIA network increases, but it is actually increasing for small values of ρ. We do not have an explanation for this observation. The execution times for the flow-level simulation are not reported in details. It varies between one second and several minutes on a machine with characteristics: processor Intel(R) Core(TM)2 Duo of 2.00GHz, 2GB RAM running Fedora Core 5 (slightly less power than the machine used for packet-level simulations). The simulation time per sample increases with the load of the system, due to the fact that the load sharing must be re-computed every time a flow starts or stops. For a given load, it also increases with the number of nodes.

We conclude that the flow-level algorithm is very efficient in terms of time. However, one question remains: how good is it in term of accuracy?

To answer this question, we first depict in Figures 1-10 the empirical complementary cumulative distribution function (CCDF) of the block download time obtained form both simulators, and for all the experiments. We report then in Table 2 the expected block service time obtained from both simulation levels and from the PS formula [START_REF] Dandoush | A realistic simulation model for peer-to-peer storage systems[END_REF]. Table 2 reports as well the relative error between results, on one hand, of NS-2 and, on the other hand, of flow-level or PS models. The relative error, for instance between results from NS-2 and FLA is calculated as follow: where E[T F LA ] is the measured empirical average download time for the FLA, and E[T N S ] is the measured download time for NS, corrected by a constant value so that the minimal values for both simulators are the same. Indeed, there is a very small difference between the minimum values in all experiments shown in the figures, which is explained by the fact that the flow-level simulator ignores the delay for establishing and closing the TCP connections and propagation delays. We have not corrected this systematic error in the figures, but we have eliminated it for relative errors. The results show that for small system load, the download time predicted by the PFFLA fits exactly that of the NS-2. The relative error between the average values is very small as shown in Table 2. The average value calculated from the PS formula is also very close but the relative error between average values of PS and NS-2 is slightly larger than that between PFFLA and NS-2. This confirms that the prediction of the duration of TCP flow is accurate. Indeed, since these are long flows, the slow start phase can be easily neglected. Other phenomena which typically perturb the throughput of TCP (packet losses, buffer fluctuations) probably happen very rarely in this case. The flow sharing algorithm apparently provides a very good approximation, for average response times as well as for distributions.

RR(N S, F LA) = E[T N S ] -E[T F LA ] E[T N S ]
When ρ is relatively large, some buffers can fill up more frequently, and then some flows tend to be relatively long in the NS-2 simulation. However, the relative errors between average values of INRIA PFFLA and NS-2 are slightly more important in this case but still very small, in particular, RR is less than 2% in the symmetric case and less than 8% in the asymmetric case. It is clear from Figures 56that the distributions measured by both simulators are different, but average values turn out to be almost identical. The same observation holds for asymmetric cases, see Figures 8(b) to 10.

The accuracy of the PS approximation for the average download time is acceptable for symmetric cases up to ρ = 36%, and degrades above ρ = 50%. The accuracy for the complete distribution can be assessed on Figure 4(b) for a load of ρ = 36%. In the asymmetric cases, the approximation is bad at low loads, and very bad at large loads. The explanation for this is the following. The download of a block at a client can be slowed down by two phenomena. The first one is that a second request arrives at the node. This is taken into account by the PS model. The second one is that one TCP flow is slowed down at the server side. This requires that at least sC u /C d blocks are downloaded simultaneously from the server. In the symmetric cases, this value is 4 and the event rarely happens, even for moderate loads. In the asymmetric case, this value is 1 and this is much more frequent. See the Appendix for more comments on the PS approximation. In order to understand better the differences between flow-level and packet-level distributions, a careful analysis of the congestion avoidance mechanism in congested networks and an extension of the algorithm to account for the overloaded system (ρ around one) are the objective of ongoing research.

Another observation is that larger the network size, better the performance of PFFLA and worse the performance of PS model as illustrated by Experiments 11, 12 and 13 (resp. 14 and 15). The number of peers in these experiments are 25, 250 and 500 respectively (resp. 250 and 500) for same load and capacities. However, the relative error occurred in Experiment 14 is less than that of 13 which is, in turn, less than that of Experiment 12. Respectively, the relative error occurred in Experiment 16 is less than that of 15. Indeed, the performance does not depend only on the system load but also on the number of peers and their capacities.

Clearly, larger the buffer sizes, better the performance in real networks. To address this point, we depict in Figure 11 the CCDF of block download time for two values of the queue limit (100 and 500 packets) in NS-2 simulation for 25 peers (50 nodes), C = 1500kbps and ρ = 70%. The relative error between the two expected download time is 6.56% (159.575 seconds for 500 packets and 170.038 for 100 packets) and then the buffer size can affect the performance of the system with high load. The main problem here is that this parameter can not take large value because it is very expensive.

We conclude that PFFLA is a reliable mechanism to analyze the service response time in many non-overloaded systems based on P2P and Grid computing concepts such as Storage Systems and Grid Delivery Networks. In particular, when the size of networks is relatively large, PFFLA predictions are very accurate as long as the system is not overloaded or close to be overloaded.

Conclusion and future work

In this report, we have proposed and analyzed the PFFLA algorithm. The algorithm is quite simple and uses the concept of "Progressive-Filling" (or max-min fairness). We have implemented the PFFLA in a flow-level simulator of parallel downloads, and we have programmed a similar model A conclusion from the literature is that different RTTs do introduce some "unfairness" in bandwidth allocations. Our next step will therefore be to find a simple yet efficient modification of Algorithm 1 to handle this situation.

E(e -sV (d) ) = (1 -ρ) (s + λ) 2 e -d(s+λ) s 2 + λ(s + (s + λ)(1 -ρ))e -d(s+λ) ,
where the load factor is ρ = λd. The inversion of this Laplace transform yields the series: .

P (V (d) ≤ d + t) (7)
For every t ∈ [0, 2d], this formula reduces to:

P (V (d) ≤ d + t) (8)
= (1ρ)e

INRIA

A.2 Small load approximations

We briefly discuss now approximations that can be performed when the load or the arrival rate is small. Consider a client downloading s flows in parallel. The nominal duration of each flow is σ/s, so that the total duration is σ if the flows are not disturbed.

Assume that the arrival of flows is a Poisson process of rate λ at all nodes: at the client node and the server nodes. Given some "tagged" download request, the probability that another request interferes with it at the client is e -λ×(2d) = e -2ρ because the request interferes if it arrives less than d units of time before or after the arrival of our tagged request. The same probability holds at each server.

The result of a request interfering at the client is that the tagged download is longer. If the arrival date of the interfering request relative to the tagged request is u, the additional response time of the tagged request is d -|u|.

The result of a request interfering at the server depends on the capacity of the server. If the capacity is enough, the interference will not slow down the flow, and the response time will not change. This is the case for the symmetric cases in our experiments. If the capacity is not enough, the flow will be slowed down. Take the case where the capacity of servers is precisely that of one of the s parallel flows. This is the case for asymmetric cases in our experiments. If an interference occurs at the server, the flow will be slowed down to half its throughput. The result is then exactly the same as when two requests interfere at the client.

Suppose now that the arrival rate λ and the load ρ are small. Ignoring the events that happen with probability o(ρ), only three events are to be considered: a) no interferences; b) one single interference at the client, none at the server; c) one single interference at the server, none at the client. According to the discussion above, we can calculate the statistics of the response time T in the two situations: Large server capacity: the probability that T = d is the probability of events a) and c), since c) does not have an influence on T . This probability is: e -2ρ = 1 -2ρ + o(ρ). Minimal server capacity: the probability that T = d is the probability of event a), that is, (e -2ρ ) 2 = 1-4ρ+o(ρ). For x ∈ [0, d], P (T > d+x) = P ( b) or c) and |u| < d-x) = 4ρ(1-x/d), and E[T ] = d(1 + 2ρ). These formulas are not in accordance anymore with the P S queueing model. This explains the bad results of the approximation in Table 2 for asymmetric cases. At the same time, this suggests a possible correction for the PS approximation formula. The prediction that P (T > d) ∼ 4ρ can however be observed on Figures 7 and8

Figure 1 :

 1 Figure 1: Experiments 1 (left) and 2 (right): progressive-filling flow-level algorithm PFFLA vs Packet-level simulation NS-2

Figure 2 :

 2 Figure 2: Experiment 3: Packet-level simulation NS-2 vs progressive-filling flow-level algorithm FLA & PS for ρ = 12%, C = 1500kbps, N = 500, S B = 8MB.

  ρ = 18%, C=1500kbps, N =500, S B =8MB. ρ = 25%, C=576kbps, N =500, S B =8MB.

Figure 3 :

 3 Figure 3: Experiments 4 (left) and 5 (right): progressive-filling flow-level algorithm PFFLA vs Packet-level simulation NS-2.

Figure 4 :

 4 Figure 4: Experiments 6 (left) and 7 (right): progressive-filling flow-level algorithm PFFLA vs Packet-level simulation NS-2.

Figure 5 :

 5 Figure 5: Experiment 8: progressive-filling flow-level algorithm FLA vs Packet-level simulation NS-2 for ρ = 50%, C = 1500kbps, N = 500, S B = 8MB.

Figure 6 :

 6 Figure 6: Experiments 9 (left) and 10 (right): progressive-filling flow-level algorithm PFFLA vs Packet-level simulation NS-2.

  ρ = 12%, C d =1500kbps, Cu=384kbps, N =50, S B =8MB. ρ = 12%, C d = 1500kbps, Cu = 384kbps, N = 500, S B = 8MB.

Figure 7 :

 7 Figure 7: Experiments 11 (left) and 12 (right): progressive-filling flow-level algorithm PFFLA vs Packet-level simulation NS-2.

Figure 8 :

 8 Figure 8: Experiments 13 (left) and 14 (right): progressive-filling flow-level algorithm PFFLA vs Packet-level simulation NS-2.

Figure 9 :

 9 Figure 9: Experiment 15: progressive-filling flow-level algorithm FLA vs Packet-level simulation NS-2 for ρ = 36%, C d = 1500kbps, C u = 384kbps, N = 1000, S B = 8MB.

Figure 10 :

 10 Figure 10: Experiments 16 (left) and 17 (right): progressive-filling flow-level algorithm PFFLA vs Packet-level simulation NS-2.

Figure 11 :

 11 Figure 11: Queue size effect in Packet-level simulation NS-2 for ρ = 70%, C = 1500kbps, N = 50, S B = 8MB.

= ( 1 (- 1 ) n e -nρ 1 {t≥nd} n m=0 n m ( 2 -

 112 ρ)e -ρ ∞ n=0 ρ) m (1ρ) n-m [λ(tnd)] 2n-m (2nm)! 1 + 2λ tnd 2nm + 1 + λ 2 (tnd) 2 (2nm + 1)(2nm + 2)

  For x ∈ [0, d], P (T > d + x) = P ( b) and |u| < dx) = 2ρ(1x/d), and E[T ] = d(1 + ρ). These formulas are in accordance with Equations 8 and (6). The prediction that P (T > d) ∼ 2ρ can be observed on Figures 1 to 4. The linear behavior of the distribution of P (T > x) for x ∈ [d, 2d] is also clear on these figures.

  (a) (with ρ = 12%). The almost-linear behavior of the distribution of P (T > x) for x ∈ [d, 2d] is also clear on these figures.

Table 1

 1

				: Experiments setup
	Experiment	N /2	C d /C u	S B /S F	1/λ	ρ	samples	Required time
	number	peers	kbps	MB	sec.	%		hours
	1	25	384/384	4/1	60	6	44853	20
	2	25	576/576	8/2	39.88	12	77651	18
	3	250	1500/1500	8/2	1.536	12	9157	7
	4	250	1500/1500	8/2	1.024	18	9801	5
	5	250	576/576	8/2	1.913	25	31695	33
	6	250	1500/1500	8/2	0.734	25	37688	31
	7	250	1500/1500	8/2	0.510	36	41104	30
	8	250	1500/1500	8/2	0.367	50	41877	36
	9	250	1500/1500	8/2	0.306	60	67500	58
	10	250	1500/1500	8/2	0.262	70	279333	264
	11	25	1500/384	8/2	59.81	12	31966	23
	12	250	1500/384	8/2	5.98	12	33509	54
	13	500	1500/384	8/2	2.99	12	31817	25
	14	250	1500/384	8/2	1.99	36	57645	54
	15	500	1500/384	8/2	0.996	36	11601	6
	16	500	1500/384	8/2	0.718	50	12580	6
	17	500	2000/384	8/2	0.718	50	9000	5

Table 2 :

 2 Measurements for the PFFLA and the packet-level simulation; comparison with the PS

	model Experiment	N /2	ρ	E[T N S ]	E[T F LA ]	RR%	E[T P S ]	RR%
	number	peers	%	sec.	sec.	(N S, F LA)	sec.	(N S, P S)
	1 symm.	25	6	96.062	95.45	0.6%	95.44	0.6%
	2 symm.	25	12	135.04	136.071	-0.7%	141.59	-4.8%
	3 symm.	250	12	51.42	52.089	-1.3%	52.195	-1.5%
	4 symm.	250	18	55.465	55.96	-0.8%	56.015	-0.9%
	5 symm.	250	25	161.252	160.196	0.6%	166.132	-3%
	6 symm.	250	25	61.068	61.517	-0.7%	61.243	-0.3%
	7 symm.	250	36	73.547	73.346	0.2%	71.7692	2.4%
	8 symm.	250	50	99.501	97.75	1.7%	91.864	7.6%
	9 symm.	250	60	129.066	127.691	1%	114.83	11%
	10 symm.	250	70	176.45	180.05	-2%	153.107	13.2%
	11 asymm.	25	12	61.137	62.901	-2.8%	52.19	17%
	12 asymm.	250	12	64.738	64.935	-0.3%	52.19	19.3%
	13 asymm.	500	12	65.298	65.182	-0.18%	52.19	20%
	14 asymm.	250	36	103.70	110.231	-6.2%	71.76	30.8%
	15 asymm.	500	36	105.82	110.396	-4.3%	71.76	32.1%
	16 asymm.	500	50	143.75	149.213	-3.8%	91.865	36%
	17 asymm.	500	50	140.94	149.213	-5.9%	68.45	51.4%

RR n°7159

The question of how to handle different round trip times is left for future work.INRIA