Olivier Roussel
email: olivier.roussel@lip6.fr

Michèle Soria
email: michele.soria@lip6.fr

Boltzmann sampling of ordered structures

Boltzmann models from statistical physics combined with methods from analytic combinatorics give rise to efficient algorithms for the random generation of combinatorials objects. This paper proposes an efficient sampler which satisfies the Boltzmann model principle for ordered structures. This goal is achieved using a special operator, named box operator. Under an abstract real-arithmetic computation model, our algorithm is of linear complexity upon free generation ; and for many classical structures, of linear complexity also provided a small tolerance is allowed on the size of the object drawn. The resulting programs make it possible to generate random objects of sizes up to 10 7 on a standard machine.

Introduction

In 2004, Duchon, Flajolet, Louchard and Schaeffer [START_REF] Duchon | Boltzmann samplers for the random generation of combinatorial structures[END_REF] proposed a new model, called Boltzmann model, which leads to systematically construct samplers for random generation of objects in combinatorial classes described by specification systems. This framework has two main features: uniformity (two objects of the same size will have equal chances of being drawn) and quasi-linear time complexity.

Boltzmann samplers depend on a real parameter x and generate an object a in a given combinatorial specifiable class A with a probability essentially proportional to x |a| where |a| is the size of a. Hence they draw uniformely in the class A n of all the objects of size n in A. The size of the output is a random variable, and parameter x can be efficiently tuned for a targetted mean value. Moreover using rejection, one can obtain exact size or approximate size samplers. Efficient generation of huge objects makes it possible to adress problem of testing and benchmarking. This new approach differs from the "recursive method" introduced by Nijenhuis and Wilf [START_REF] Nijenhuis | Combinatorial algorithms[END_REF] bringing the possibility of relaxing the constraint of an exact size for the output. This implies a significant gain in complexity: no preprocessing phase is needed and expected time complexity is linear in the size of the output.

Boltzmann samplers have been developped for a whole set of combinatorial classes, labelled or unlabelled. Basically, they are classes defined from basic elements by means of fundamental constructions such as cartesian products, disjoint unions, sequences, sets and cycles, these operators being well known in combinatorics.

1

In this paper, we focus on a particular operator, introduced by Greene [START_REF] Greene | Labelled formal languages and their uses[END_REF], which allows to construct objects with internal order constraints. We extend Boltzmann model to efficiently generate these objects. The main idea is to slightly change the value of the x parameter along the execution of the algorithm. This simple idea leads to Boltzmann samplers for objects such as alternating permutations (which we could not handle before) and more generally any variety of increasing tree [START_REF] Bergeron | Varieties of increasing trees[END_REF]. This also leads to a new point of view for random sampling of combinatiorial objects built upon Cycle, Set or Sequence, as these constructions can be rephrased using ordered structures [START_REF] Greene | Labelled formal languages and their uses[END_REF][START_REF] Flajolet | Analytic combinatorics[END_REF]. Nonetheless, there is still some points to be adressed in order to make this method fully effective: we have to evaluate the generating series of our objects in some arbitrary value, these generating functions may not exist in a closed form, and we may have to solve differential equations.

As an example, we will focus on alternating permutations. All the proofs are available in the appendix.

Boltzmann model for labelled structures

Definition A combinatorial class A is a denombrable (or finite) set, with a size function size | • | : A → N and such that there is only finitely many objects of each size.

In the following, we will use these notations: if A is a class, then for any object α ∈ A, |α| is its size. Furthermore

A n = {α ∈ A | |α| = n} and a n = Card(A n ).
There exists two kind of combinatorial objects: unlabelled and labelled ones. Basically, a labelled object of size n could be seen as an unlabelled object of the size, in which each atomic part is numbered by a number in {1, . . . , n}. In other words, a labelled object of size n is simply an unlabelled one of the same size paired with a permutation in S n .

In all the following, we will only consider labelled objects.

Definition Let A be a (labelled) combinatorial class. We define the exponential generating function defined by

A(z) = α∈A 1 |α|! z |α| = n∈N a n z n n!
Now, in order to construct some objects, we need a set of rules. Precisely, we will have some basics objects, called atoms, and a set on operators which allow us to build big objects from littler ones. A few of these operators are presented on figure 1 or more completely in [START_REF] Flajolet | Analytic combinatorics[END_REF].

Definition A Boltzmann sampler ΓC(x) for a (labelled) combinatorial class C is a random generator such that the probability of drawing a given object γ ∈ C of size n is exactly:

P x (γ) = 1 C(x) x |γ| |γ|! = 1 C(x) x n n! where C(x) is obviously the (exponential) generating function of C. A Description A(z) ΓA(x) ε Empty class 1 return ε Z Atomic class z return ➀ B × C Cartesian product B(z) × C(z) return (ΓB(x), ΓC(x)) B + C Disjoint union B(z) + C(z) if Bernouilli B(x) B(x)+C(x)
then return ΓB(x) else return ΓC(x) Seq(B) Sequence Let's notice that this sampler depends on a (real) parameter x, which can be tuned to aim a given expected size for the output of the algorithm. To be precise, if we denote by N the size of the output, one can solve the équation

E x (N ) = x C ′ (x)
C(x) in order to fix the value of x. We could also compute then

V x (N ) = x 2 C ′′ (x) C(x) + x C ′ (x) C(x) 1 -x C ′ (x) C(x)
. Boltzmann model is a very useful tool to efficiently generate combinatorial structures [START_REF] Duchon | Boltzmann samplers for the random generation of combinatorial structures[END_REF][START_REF] Flajolet | Boltzmann sampling of unlabelled structures[END_REF]. In particular, it is possible to automatically build a sampler according to the specification of a combinatorial class, following recursively the rules described in figure 1 or in [START_REF] Duchon | Boltzmann samplers for the random generation of combinatorial structures[END_REF].

Extension to ordered structures

We introduce the central object of this paper: the box operator, for which we extend the Bolzmann formalism.

The main idea of this paper consist of changing the value of the x parameter of the sampler ΓA during the recursive calls. Indeed, one can notice that, in the classical Boltzmann sampler (see figure 1), this value is fixed at the beginning and stays constant during the execution of the algorithm. Here, we will change this value during the execution according to a given probability density.

Definition The box product of two labelled classes B and C, noted by B × C, is the subset of B × C such that the least label is in the left composant of the pair, on the B part.

Furthermore, we need to ensure that B 0 = 0, i.e. the class B has no empty objects.

We have to highlight the fact that this operator is a binary one, and not an unary one! Indeed, the real operation is not • , but rather • × •.

Proposition 2.1. The generating function of

A = B × C is A(z) = z 0 B ′ (t)C(t)dt

Boltzmann sampler for the box operator

As said before, we want to change the value of the parameter x of the sampler according to a given probability law.

Proposition 2.2. If 0 < x < ρ A , then δ A x (t) = xB ′ (tx)C(tx) A(x)
is a probability density over [0, 1]. Now, we can almost explicitely present our algorithm for generating A = B × C following Boltzmann's ideas. We will need just one last notion: the derivative class.

Definition

The derivative class of a given combinatorial class B, written B ′ , is the set of objects of B in which an atom has been substituted by a hole, by a reservation. In other words, one atom (of size 1 by definition) has been replaced by a distingued atom, written ⊙, whose size is 0.

We will explain later how to compute this differentiated class. The proof just consists in checking that we can get all the desired objects, only these ones, and each one with the correct probability.

Nonetheless, we still have several points to adress before being able to really implement this algorithm: we have to define the differentiated class B ′ , and explain how to draw according to the probability density δ x .

Theoretical complexity

The free generation of objects (with no constraint on the size) works according to the following algorihm:

1. First, draw a random number according to δ x 2. Then generate two sub-objects 3. Eventually re-label the final object At this point, we first investigate the complexity of this algorithm related to the size of the generated object. Let us assume that the first step is in O(1). We will see later that this assumption is indeed correct, even if the hidden constant can be quite big. Moreover, we will suppose that the third step is also in O(1).

Then, according to the algorithm, to get an object of size n, we have to get two sub-objects of sizes k and n -k. From this, we check that the complexity of our algorithm is O(n): it's the same complexity as for the previously known Boltzmann samplers.

In a majority of real-world usage of random sampling, we just care about the shape, the form of the squeleton of the huge objects we get. The exact position of the label is not as important as the structure of the objects we get.

In addition to this complexity of the free generator, one can estimate the complexity for approximate-size generator. First, let us remind that approximatesize generation means that we want to get an object of size in

[(1 -ε)n, (1 + ε)n]
In the original paper [START_REF] Duchon | Boltzmann samplers for the random generation of combinatorial structures[END_REF], it is proved that, under some conditions on the analytical nature of the generating function, we only need a constant number of trials for approximate-size sampling. Here, these conditions are often true, and the generation stays linear.

Proposition 2.4. For any combinatorial structure whose generating function ∼ C(1 -x) -α log β 1 1-x , using a simple rejection method, the number of trials is constant for a given precision.

As a majority of "interesting" classes ave thus property, it means that we can quite often ensure a linear complexity even for approximate-size generation.

For example, in our particular case of alternating permutations, if ε = 10%, then the average number of trials is less than 7.

Example: Alternating permutations

Alternating permutations (also known as increasing proper binary trees) are defined by the specification T = Z + Z × (T × T ). From this one we can get

T (z) = z + z 0 dz dz (u)T 2 (u)du, and therefore T (z) = tan(z) with ρ T = π 2 . To lighten, let us note A = Z × T 2 and A(x) = tan(x) -x.
From there, one can write:

δ A x (t) = x dz dz (tx)T 2 (tx) A(x) = x tan 2 (tx) tan(x) -x
Hence, a Boltzmann sampler ΓT (x), for 0 x < ρ T = π 2 , is the following:

Algorithm 2 Boltzmann sampler for T Input: x ∈ R Output: An object of T , meaning an alternating permutation Require:

0 x π 2 1: Draw uniformly U over [0, 1] 2: if U x T (x) = x tan(x) then 3:
T ← Leaf ➀ 4: else 5:

Draw α over [0, 1] according to δ A x 6:

T ← Node ⊙ (ΓT (αx), ΓT (αx)) 7: end if 8: return T with right labels One can compute the expected size of the output:

E T x (N ) = x T ′ (x) T (x) = x 1 + tan 2 (x) tan(x) = 2x sin(2x) ∼ π/2 1 1 -x π/2
More precisely, one can compute the distribution of the sizes of the output. We have P T

x (N = 2n) = 0 and a simple computation leads to

P T x (N = 2n + 1) = T 2n+1 x 2n+1 T (x)(2n + 1)! ∼ ∞ 4 π tan x x π/2 2n+1
This probability asymptotically follows a geometric distribution of parameter x π/2 . This behaviour really differs with the proper binary tree case.

Algorithmic issues

In order to be able to implement on a computer the previous algorithm, we first have to adress a few practical points.

The B ′ combinatorial class

This class B ′ intuitively corresponds with the objects of the B class, in which a reservation replaced an atom. One can check that the generating function of

B ′ is (B ′ ) (z) = dB(z) dz (z).
Proper definition and set of properties are presented in [START_REF] Bergeron | Combinatorial species and tree-like structures[END_REF].

To define it precisely, we have to introduce a new special atom, written ⊙ and such that | ⊙ | = 0. But be aware that ⊙ = ε ! This atom, simply speaking, is juste a reservation for a future atom Z. It's just an empty place, but an existing and reserved one.

Then, one can recursively define or class B ′ by:

ε ′ = ∅ ⊙ ′ = ∅ Z ′ = ⊙ (A + B) ′ = A ′ + B ′ (A × B) ′ = A ′ × B + A × B ′ (Seq(A)) ′ = Seq(A) × A ′ × Seq(A) (Set(A)) ′ = A ′ × Set(A) (Cyc(A)) ′ = A ′ × Seq(A) (A × B) ′ = A ′ × B
We can check that we can always compute a formula without any differentiated class isomorph to our differentiated class. Hence, one can automaticaly compute a Boltzmann sampler for this isomorph class, and so for the differentiated one. So, we are able to draw in a linear time an object of this class, following Boltzmann's model.

This being said, we have to notice that the classes B and B ′ are quite close. This observation beeing done, it is really frustrating to have to unroll this whole rewriting system to be able to draw objects in the differentiated class. It might be possible to find a more efficient and satisfying way to draw into such classes.

Drawing according to δ x

The second issue is that we have to be able to draw a (real) variabel according to an (almost arbitrary) distribution δ x over [0 , 1] One can easily verify that, if x < ρ, then δ x is a C ∞ function from [0, 1] to R + . More precisely, we know that ∀k ∈ N, ∀t ∈ [0, 1], ∂ To draw a random variable X according to δ x , we will adapt a rejection method to our particular distribution. We define a subdivision of [0, 1] into consecutive intervals ([α 0 , α 1 ], . . . ,

[α i-1 , α i ], . . . , [α N -1 , α N ]) with 0 = α 0 < α 1 < • • • < α i < • • • < α N = 1.
On each [α i-1 , α i ], let h i δ x be a continuous function. We choose the h i simple enough to be able to compute and know all the values A i = (α i -α i-1 ) αi αi-1 h i (u)du. These conditions on the h i ensure that we will be able to draw a random variable according to these functions, and to correct the fault induced by this approximation of δ x .

In the following, let Area(h) = N i=1 A i be the total area under all the h i . We use then the following algorithm: Algorithm 3 Draw a variable x according to a probability density δ x Input:

δ x : [0, 1] → R + ; N ∈ N ; α i ∈ [0, 1] ; and h i : [α i-1 , α i ] → R + Output: X ∈ [0, 1]
Require: δ x is a probability density over [0, 1] ; α 0 = 0 and α N = 1 ; and ∀i ∈ {1, . . . , N }, h i δ x Ensure: X ∼ δ x 1: repeat 2:

Draw i ∈ {1, . . . , N } such that P(i = k) = A k Area(h) 3:
Draw a random variable X according to the density hi Ai sur

[α i-1 , α i ] 4:
Draw a random variable Y according to the uniform law over [0, 1] 5: until Y h i (X) δ x (X) 6: return X Theorem 3.1. This algorithm returns a random variable X according to the distribution δ x ; and the expected number of rounds to finish is Area(h). We can compute that, if p = 1 Area(h) , then P(stop k ) = (1 -p) k-1 p and

0 1 t δ A x (t)
E({Number of trials}) = ∞ k=0 kP(stop k ) = 1 p = Area(h).
We can see here that the h i functions have to be the closest to the density δ x , in order to minimize the average number of trials to draw the variable X.

We tried two kind of approximations for the h i functions:

• piecewise constant functions:

∀i ∈ {1, . . . , N }, ∀t ∈ [α i-1 , α i ], h i (t) = δ x (α i )
• piecewise linear functions δ x being non-decrasing, we are sure that in both case:

∀i ∈ {1, . . . , N }, ∀t ∈ [α i-1 , α i ], h i (t) δ x (t)
In practice, we can see that for our problem, piecewise constant functions appear to offer the better compromise between pre-computing time and computing time. As our distributions are "more and more increasing" over [0, 1]i.e. all their derivatives are non-negatives -we chose to have a geometric subdivision of the interval: α i = 1 -2 i and α N = 1 (see figure 2). Experimentally, this kind of subdivision works quite well.

Conclusion

The random generation of ordered structures is in general a very difficult problem. In this short paper, we have presented a way to efficiently generate some of these structures, using Boltzmann model. As said before, it is possible with our algorithm for the box operator to generate objects of size up to 10 7 in reasonable time (at most about 100 seconds).

Moreover, we keep all the good properties induced by Boltzmann samplers: genericity of the algorithms, uniformity over objects of the same size, linear complexity in the size of the output as far as only the shape is concerned. The basic idea behind this extension is to slightly change the parameter during the recursive calls.

We think this idea can be extended to generate a lot of other structures and operators. For example, the shuffle operator, the Hadamard product, or almost any specification which can be described by a system of differential equations on the combinatorial classes.

Proof. By definition of δ A

x , we can write that ∀t ∈ [0, 1], 0 δ A x (t), and:

1 0 δ A x (t)dt = 1 0 xB ′ (tx)C(tx) x 0 B ′ (u)C(u)du dt = x x 0 B ′ (z)C(z) dz x x 0 B ′ (u)C(u)du = 1 Therefore, δ A
x define a correct probability density over [0, 1].

Theorem The algorithm described above for the combinatorial class A = B × C is a correct Boltzmann sampler.

Proof. We have to check that we can get each object in A with this algorithm, and only those ones ; and mainly that we get each object with the right probability.

The first point is easy to check, as the algorithm follows by construction the exact definition of what is a box product.

About the second point, we have to check that we get any object α ∈ A with probability P

A x (α) = x |α| A(x) 1
|α|! . Indeed, we can remind that:

C(x) = γ∈C x |γ| |γ|! and B(x) = β∈B x |β| |β|!
From this equation:

B ′ (x) = β∈B x |β|-1 (|β| -1)! A(x) = x 0 B ′ (u)C(u)du = x 0   β∈B u |β|-1 (|β| -1)!     γ∈C u |γ| |γ|!   du = x 0 (β,γ)∈B×C u |β|+|γ|-1 (|β| -1)!|γ|! du = x 0 α=( e β,e γ)∈A (| β| -1)!| γ|! (|α| -1)! u |α|-1 (| β| -1)!| γ|! du = x 0 α∈A u |α|-1 (|α| -1)! du A(x) = α∈A x |α| |α|!
Then we can write, assuming having Boltzmann samplers for B ′ and C, for a given pair (β ′ , γ) ∈ B ′ × C, the probability P

x ((β ′ , γ)) = x |β ′ | B ′ (x) 1 |β ′ |! x |γ| C(x) 1 |γ|! .
Hence, for a given α = ( β, γ) ∈ A, we can write, for the algorithme described earlier for A that:

P A x (α) = 1 0 (| β| -1)!| γ|! (|α| -1)! P tx ((β ′ , γ)) δ A x (t)dt = 1 0 (|β| -1)!|γ|! (|α| -1)! × (tx) |β ′ | B ′ (tx) 1 |β ′ |! (tx) |γ| C(tx) 1 |γ|! × xB ′ (tx)C(tx) x 0 B ′ (u)C(u)du dt = x |α| A(x) 1 (|α| -1)! 1 0 t (|α|-1) dt = x |α| A(x) 1 |α|!
Eventually, we get a correct Boltzmann sampler for this combinatorial class.

Theorem The algorithm 3 returns a random variable X according to the distribution δ x ; and the expected number of rounds to finish is Area(h).

Proof. It is quite simple to see that this algorithm implements the rejection method for the distribution δ x under the curve defined by the h i , and since h i δ x on each inverval, the algorithm is correct.

We will now find the expected number of rounds before returning a value. Let us use the events {The algorithm ends after exactly k rounds} written stop k , and {We chose the i th interval} written A i . First, we will have to compute P(stop 1 ) = N i=1 P(stop 1 |A i )P(A i ) But P(A i ) = Ai Area(h) according to the algorithm. And as to the conditional probability, we get:

P(stop 1 |A i ) = P(Y h i (X) δ x (X)|A i ) = αi αi-1 P(Y h i (X) δ x (X)|X ∈ [µ, µ + dµ], A i )P(X ∈ [µ, µ + dµ]|A i ) = αi αi-1 P(Y h i (X) δ x (X)|X ∈ [µ, µ + dµ], A i ) h i (µ) A i dµ = αi αi-1 δ x (µ) h i (µ) h i (µ) A i dµ = 1 A i αi αi-1

Alternating permutations: experimental results

On a 1Ghz notebook, with 2Go of RAM, and a naive implementation in OCaml, we can reach sizes about 10 7 in reasonable time (about 10 seconds for free generation: see figure 4). Moreover, on can compute, for this specific class, the number of trials for approximate size generation. We can show that, for a given precision ε, only Using this approach, we can use Boltzmann model advantages, on so draw large object in very few time. For example, the example displayed as figure 3 can be computed in less than one millisecond on a standard computer.

1 1 -

 1 B(z) l := Geom(B(x)) return (ΓB(x), . . . , ΓB(x)) l times Set(B) Set e B(z) l := Poisson(B(x)) return {ΓB(x), . . . , ΓB(x)} l times Cyc(B) Cycle log 1 1-B(z) l := Logarithmic(B(x)) return ΓB(x), . . . , ΓB(x) l times

Figure 1 :

 1 Figure 1: Some classical constructors, their generating function and sampler

Algorithm 1 2 : 3 : 4 :

 1234 Boltzmann sampler ΓA for A = B × C Input: One real number x Output: An object of A Require: 0 < x < ρ A 1: Draw U in [0, 1] according to the probability density δ A x Randomly draw an object b ′ from the class B ′ using ΓB ′ (U x) Randomly draw an object c from the class C using ΓC(U x) Let a = (b ′ , c) correctly labelled 5: return The object a, with the reservation ⊙ replaced by an atom Z One should notice that the shape of this algorithm really looks like the one for the simple cartesian product B × C. Particularly, both of the two recursive calls for B ′ and C are independant and both with the same modified value for the parameter.

Theorem 2 . 3 .

 23 The algorithm described above for the combinatorial class A = B × C is a correct Boltzmann sampler.

  k δx ∂t k (t) 0. Particularly, the probability density is non-decreasing on his definition domain, and ∀t ∈ [0, 1], lim x→0 δ x (t) = 1 and lim x→ρ δ x (t) = Dirac(t -ρ).

Figure 2 :

 2 Figure 2: Shape of a typical δ x distribution

Figure 3 :

 3 Figure 3: Increasing proper binary tree: size 2033, drawn in less than 1ms

  From this point, P(stop k ) = (1 -p) k-1 p and E({Number of rounds}) =

e 4 1 sinhFigure 4 :

 14 Figure 4: Time for generating alternating permutation

A(x)

Proof. Following the definition of the box product, we have

Indeed, we have to distribute only n-1 labels, the least being fixed in B. Among those, we have to choose only k -1 for the first composant, and n -k for the second one. We can then rewrite this equation, using properties of binomial coefficients, and the fact that B 0 = 0:

which is a Cauchy product. This last remark leads to the claimed formula.

Remark Funnily, using this framework, the well-known integration by parts formula

can be proven and read as: "In a pair of labelled object, the least labeal is either on the first composant, either on the second one", which is absolutely trivial.

Proposition If 0 < x < ρ A , then δ A x (t) = xB ′ (tx)C(tx)


PhD Nicolas Cedilnik
email: nicolas.cedilnik@inheart.fr

PhD Jean-Marc Peyrat

Weighted tissue thickness

Keywords: imaging, thickness, infarct, fat

Measuring the thickness of a tissue can provide valuable clinical information; anatomical structures segmented on medical images can include sub-structures (inclusions) corresponding to a dierent biological tissue. This article presents a method, based on partial dierential equations, to measure the thickness of one specic tissue in this particular conguration.

After describing the mathematical formulation of our weighted thickness denition, we show on synthetic geometries in one, and two dimensions that it outputs the expected results. We then present three possible applications of our method on cardiac imaging data: measuring the muscular thickness of a ventricle with fat inltration; measuring the thickness of an infarct scar; visualising the transmural extent of an infarct scar.

Introduction

Measuring the thickness of an anatomical structure is a common task for radiologists. More specically, in the eld of cardiology, it is often needed to measure the thickness of the myocardial muscle, as it convey information about its health status, which has been shown to be related to its electrophysiological properties ( [START_REF] Michael | Correlation between computer tomography-derived scar topography and critical ablation sites in postinfarction ventricular tachycardia[END_REF], [START_REF] Takigawa | Detailed comparison between the wall thickness and voltages in chronic myocardial infarction[END_REF]).

The visible myocardium segmented on medical images can have inclusions of non-muscle tissue such as fat inltration, calcications, or brosis. It can happen that we are interested in measuring the thickness of the actual muscle bers, ignoring these inclusions (g. [START_REF] Cedilnik | VT scan: Towards an ecient pipeline from computed tomography images to ventricular tachycardia ablation[END_REF]. In this article we describe a method to perform this task, and show that it can reciproqually be used to measure the thickness of the inclusion.

Methodology

Our method is an extension of the method described by Yezzi and Prince ( [START_REF] Yezzi | An eulerian PDE approach for computing tissue thickness[END_REF]).

R: myocardium blood pool inclusions

∂ R ∂ R 0 1
th ic k n e s s Fig. 1. Illustration of the problem adressed by our method. We want to measure the thickness of R, excluding the thickness of the inclusions, ie the length of the thickness line, exclusing the dotted part. ∂0R: inner boundary. ∂1R: outer boundary.

Correspondence trajectories

Briey, Yezzi and Prince dene thickness at each point x of the tissue region R as the total arclength of a unique curve, passing through x. This curve originates on the the inner boundary of the tissue region ∂ 0 R, and terminates on its outer boundary ∂ 1 R (see g. 1)

Such curves can be obtained by solving the Laplace equation over R:

∆u = 0 (1)
with the Dirichlet boundary conditions:

u(∂ 0 R) = 0 and u(∂ 1 R) = 1 (2)
u denes a scalar eld in R and is used to dene curves within R with the tangent eld

- → T : - → T = ∇u ||∇u|| (3)
For didactic purposes, in g. 1 the curve is represented as a straight dotted line (thickness), but -→ T actually denes curvy trajectories for complex shapes of R.

Arclength computation

Yezzi and Prince ( [START_REF] Yezzi | An eulerian PDE approach for computing tissue thickness[END_REF]) dene the length functions L 0 , where L 0 (x) gives the arclength of the correspondence trajectory between ∂ 0 R and x (reciprocally, L 1 (x) between ∂ 1 R and x).

∇L 0 • - → T = 1, with L 0 (∂ 0 R) = 0 (4)
Weighted tissue thickness 3

-∇L 1 • - → T = 1, with L 1 (∂ 1 R) = 0 (5)
Thickness at every point x is then dened by summing these length functions:

W (x) = L 0 (x) + L 1 (x) (6)

Weighted arclength

To account for the inclusion of tissue which thickness we want to exclude, we propose to replace the right term (constant) of eq. ( 4) and eq. ( 5) by a function f varying over the tissue region.

∇L 0 • - → T = f (x) (7) -∇L 1 • - → T = f (x) (8)
We call f the arclength weight function and it has the following properties:

If x is fully occupied by the tissue region we want to measure, f (x) = 1 If x is fully occupied by an inclusion we want to exclude from the thickness measure, f (x) = 0 f (x) can also take any value in the [0, 1] interval. This can be useful to account for partial volume eects by using a transfer function between the intensity (eg, Hounseld Units) of the medical image and f .

Ray-tracing as a possible alternative

Another approach to achieve comparable measurements is to use ray-tracing. In such framework, one needs to dene straight lines emanating from one surface and measure the length of the segment between the two surfaces ∂ 0 R and ∂ 1 R.

Typically, these lines would be normal to the surface of the considered anatomical structure at the point where they emanate. The measured distance (ie, the tissue thickness) could be similarly weighted with the weight function f . However, our methodology presents several advantages over this approach:

It operates directly on voxels and it does not require to dene a surface, something usually achieved by converting voxel data to triangular meshes.

While there are several algorithms and tools to achieve this conversion, arbitrary choices must be made to control, among other things, the smoothness of the surface; and these choices, inuencing the directions of the rays, can have drastic eects on the thickness values.

The results with our method are thickness values for all voxels of the segmented structure, and not on a single surface. This can be leveraged for further analysis, or parameterizing fast electrophysiological models based on rectangular-grid data ([4], [START_REF] Cedilnik | VT scan: Towards an ecient pipeline from computed tomography images to ventricular tachycardia ablation[END_REF]).

Nicolas Cedilnik, PhD and Jean-Marc Peyrat, PhD

There is no need to dene out of which surface these trajectories have to be calculated. Since normals have dierent directions depending on whether they are dened from the inner or outer surface, this would result in visually dierent thickness maps depending on the considered surface. With our method, and as a consequence of both the denition of the correspondence trajectories and the voxel nature of the output, one can see topologically similar thicknesses on the epicardial, endocardial, or even mid-wall surface.

3 Implementation

Algorithm

It has been shown ( [START_REF] Yezzi | An eulerian PDE approach for computing tissue thickness[END_REF]) that computing eq. ( 4) and eq. ( 5) over R, in 3D amounts to solving (equations ( 8) and ( 9) of [START_REF] Yezzi | An eulerian PDE approach for computing tissue thickness[END_REF])

L 0 [i, j, k] = f (x) + |T x | L 0 [i ∓ 1, j, k] + |T y | L 0 [i, j, ∓1, k] + |T z | L 0 [i, j, k ∓ 1] |T x | + |T y | + |T z | L 1 [i, j, k] = f (x) + |T x | L 1 [i ± 1, j, k] + |T y | L 1 [i, j, ±1, k] + |T z | L 1 [i, j, k ± 1] |T x | + |T y | + |T z | NB: in [6], f (x) is a constant equal to 1.

Numerical solving

Results can be obtained by the same methods described by Yezzi et Prince ( [START_REF] Yezzi | An eulerian PDE approach for computing tissue thickness[END_REF]). Iterative approaches until convergence are possible, but a fast-marchinglike algorithm is more ecient (computational resource-wise), especially in a single-threaded context, since it requires a single pass over R, in theory. However we want to note here, than unlike [START_REF] Yezzi | An eulerian PDE approach for computing tissue thickness[END_REF] describe, our experiments showed that this fast-marching-like approach does not reach convergence in a single pass. We propose an iterative/fast-marchine-like combined approach where the traversal order of the rst pass is memorized for subsequent passes, reaching convergence more rapidly than a naive iterative approach. 1

Results on toy geometries

One dimension We conducted experiments in one dimension to verify that the results were those expected. We set correspondence trajectories to straight horizontal lines and set a 4-element wide hole along these trajectories where the weight function values ranged from 0 to 1 (g. 2).

1 An informal benchmark with a real left ventricular wall geometry of 363268 voxels (voxel spacing = 0.8 mm), on an AMD Ryzen 9 3950X CPU, single-threaded, took 11.2 ± 0.2 seconds for the iterative approach vs 7.5 ± 0.2 seconds for the semi-ordered approach over 10 runs. 3)); weight function f ; weighted arclengths L0 and L1 (eq. ( 4), eq. ( 5))

Diffusion

and the resulting weighted thickness.

As expected, the resulting weighted thickness W p (x) equals to the number of elements n in R along the horizontal line when the weight function f (x) is set to 1 along the line, corresponding to the non-weighted thickness computation.

When f (x holes ) = 0, W p (x) = n -n holes . For 0 < f (x holes ) < 1, n < W p (x) < -n holes .

Diffusion

Weight L0 L1 Thickness

Applications

In this section we present possible applications of our method on clinical data.

Workow

For all clinical applications, the preliminary step is a segmentation of the structure that we want to measure. For the myocardium, this is typically done by segmenting the blood pool inside the cardiac cavities to create an endocardial binary mask, and the outer layer of the myocardial wall to create an epicardial mask. [START_REF] Komatsu | Regional myocardial wall thinning at multidetector computed tomography correlates to arrhythmogenic substrate in postinfarction ventricular tachycardia: Assessment of structural and electrical substrate[END_REF] One then has to dene the weight function f . The way it is dened varies depending on the application, as shown in the following sections. Our method can also be used to quantify the local thickness of an infarct scar.

To achieve this, after segmentation of the endocardial and epicardial masks, one has to segment the scar mask (how to obtain such segmentation is outside the scope of this article).

The weight function f is then dened such that f (x) = 0 in the myocardium and f (x) = 1 in the scar. This can be used to visually assess the thickness of the brosis, as shown on g. 5

Scar transmurality

Instead of the absolute thickness of the scar, in millimeters, our method can also be used to visualise the local proportion of brosis in the myocardial wall, a ratio we call scar transmurality. A value of 1 would mean that an area of the myocardial wall is pure brosis and reciproqually, a value of 0 must be interpreted as the absence of brosis in this area, ie. no ischemic scar.

To obtain such measure, a possible approach is to:

1. Compute the non-weighted thickness of the wall, ie, setting f (x) = 1 all over the myocardial wall. The transmural extent, i.e., the scar transmurality is dened for every voxel of the myocardial wall. See g. 6 for an example result.

Conclusion

In this article, we dened a measure we named weighted thickness, suited to measure the thickness of a tissue on medical images, when it presents inclusions of a dierent tissue that we want to exclude from the measure. It does so by giving voxels dierent weights depending on a weight function that must be dened accordingly to the considered application. Like other partial dierential equationbased thickness measures, it outputs a scalar map of the same dimensions as the input segmentation, in which each pixel (or voxel) has a thickness value, coherent from the inner to the outer boundary of the measured structure.

The main limitation of our article is that we do not evaluate the clinical relevance of the weighted thickness measure. In particular, the transfer function described in 4.2 has been arbitratily dened and would require tuning and validation for the implied application, ie dening the health status of the myocardial on CT images. Our future work will focus on such validation.

Fig. 2 .

 2 Fig.2. Toy example in 1D. From top to bottom, each row has an inclusion with a increasing weight. From left to right: scalar eld u (diusion) from eq. (1) and tangent eld -→ T (eq. (3)); weight function f ; weighted arclengths L0 and L1 (eq. (4), eq. (5))

Fig. 3 .

 3 Fig. 3. Toy example in 2D. Refer to the legend of g. 2 for details.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Myocardium-only thickness. From left to right: (top) CT image showing fat inclusions in the left myocardium (red arrows); (bottom) result of f (x); non-weighted thickness; weighted thickness showing more severe thinning. The thicknesses are shown both on a 2D slice (top) and projected on a surfacic mesh of the ventricular wall.

Fig. 6 .

 6 Fig. 6. Scar transmurality. From left to right, total, total wall thickness, scar thickness, scar transmurality. Top images are overlay of the thickness values over the original MR image, bottom images are projections of the values on a surfacic mesh.


Tegawendé Bissyandé
email: bissyand@enseirb.fr

Laurent Réveillère
email: reveillere@labri.fr

David Bromberg
email: bromberg@labri.fr

UbiGate: A Gateway to Transform Discovery Information into Presence Information

Keywords: UbiGate, SIP, Presence, Pervasive, SDP, Service Discovery

Pervasive computing involves various entities which need to coordinate tasks and share resources through different service discovery protocols. However, the multiplicity and the incompatibility of those protocols have made interconnectivity problematic. Moreover, most service discovery protocols require a strong participation of users to genuinely play their part. Consequently, service discovery in a pervasive environment has become a challenge that researchers as well as practicionners have tried to overcome through various approaches. Nevertheless, existing solutions mostly consist of designing new protocols which usually address specific application needs while participating in the increase of heterogeneity.

To address these problems, we present a new paradigm for service discovery involving the use of a gateway, called Ubi-Gate, and relying on SIP, a widespread signaling protocol. Centered around the notion of presence, UbiGate enables real time availability of service information while hiding the heterogeneity of underlying protocols. We have developped a prototype of UbiGate supporting service discovery protocols such as the protocol used in Bluetooth service discovery mechanism and a protocol enabling the detection mechanism of RFID. Preliminary results show that UbiGate enables new service discovery protocols, either IP or non-IP based, to be seamlessly supported with no significant overhead in discovery latency.

INTRODUCTION

Ubiquitous Computing, otherwise known as pervasive computing, refers to a computing paradigm where the user is constantly in interaction with computing devices [START_REF] Kindberg | System software for ubiquitous computing[END_REF]. As envisioned by Weiser, pervasive computing takes into account the natural human environment and allows the computers themselves to vanish into the background [START_REF] Weiser | The computer for the 21st century[END_REF]. In such environments, computers surrounding the user share communication resources and computing capacities packaged as services. The pervasive computing scheme therefore aims at enabling the use of those services in the most seamless and effortless possible way. Users collect information that describes attributes of available services using Service Discovery Protocols (SDPs). Currently, pervasive computing makes use of a wide range of SDPs [START_REF]Specification of the Bluetooth System Core version 1.0b: Part E, Service Discovery Protocol (SDP)[END_REF][START_REF] Guttman | Service Location Protocol, Version 2[END_REF][START_REF] Helal | Konark -a service discovery and delivery protocol for ad-hoc networks[END_REF][START_REF] Liang | Mobile service discovery protocol (msdp) for mobile ad-hoc networks[END_REF][START_REF] Marin-Perianu | Cluster-based service discovery for heterogeneous wireless sensor networks[END_REF][START_REF] Mian | A survey of service discovery protocols in multihop mobile ad hoc networks[END_REF][START_REF]Jini Architecture Specification version 2.0[END_REF][START_REF]Salutation Architecture Specification version 2.0c[END_REF][START_REF]UPnP Device Architecture version 1.0[END_REF]. Therefore, a user will be isolated if his device lacks the capacity to discover resources available in its vicinity because of protocol incompatibility.

Existing SDPs are either based on the Pull model or the Push model. In the Pull model, service requesters regularly send requests to discover services available in their environment. In the Push model, service providers automatically announce the services they provide by sending advertisements. However, neither model correctly address the constraints of pervasive computing. Indeed, the volatility of devices requires service requesters to reiterate their discovery requests to get up-to-date information when using the Pull Model. In the Push Model, advertisements sent by service providers are performed cyclicly. Hence, users must wait for the next advertisement to get accurate information. The dynamicity of pervasive environments suggests that users should be notified when the state of a service changes in order to provide realtime snapshots of the context in which they operate. Such a realtime service discovery can be compared to the notification of contact status in Instant Messaging (IM) clients such as Windows Live T M Messenger or Google Talk T M . In those IM clients, the user performs a realtime monitoring on the presence status of each member of its buddy list. This ability to access realtime information about a person's status, communications capabilities, and preferences as suggested by Rosenberg in [START_REF] Rosenberg | Presence: The best thing that ever happened to voice[END_REF] is relevant to virtually every means of communication.

So as to convey presence information dynamically, various instant messaging protocols are available, such as XMPP (also known as Jabber) [START_REF] Saint-Andre | Extensible Messaging and Presence Protocol (XMPP): Core. RFC 3920[END_REF], IRC [START_REF] Oikarinen | Internet Relay Chat Protocol[END_REF] and SIP/SIMPLE [START_REF] Lonnfors | Session initiation protocol (SIP) extension for partial notification of presence information[END_REF][START_REF] Roach | Session Initiation Protocol (SIP)-Specific Event Notification[END_REF][START_REF] Rosenberg | A presence event package for the session initiation protocol (SIP). RFC 3856[END_REF]. In particular, the Session Initiation Protocol, otherwise known as SIP, is widely adopted for its simplicity, portability and extensibility [START_REF] Jiang | Towards junking the PBX: Deploying IP telephony[END_REF][START_REF] Rosenberg | Presence: The best thing that ever happened to voice[END_REF].

Taking into consideration both the aforementioned challenges of Service Discovery in pervasive computing and the SIP widespread adoption, we introduce a new service discovery paradigm based on SIP.

This paper

This paper presents a new approach for service discovery based on the SIP protocol. It consists of exposing discovery information about services as presence information for users. To do so, we introduce the UbiGate gateway which is in charge of transforming discovering information into presence information. UbiGate handles all discovery tasks in the environment and delivers pertinent information to each user.

The contributions of this paper are as follows:

• We have proposed a new approach for service discovery involving the use of a gateway and relying on the widely used SIP protocol. Following our approach, we present UbiGate: a gateway to transform discovery information into presence information;

• We have developped dedicated modules for UbiGate to support RFID and Bluetooth technologies. Our preliminary experiments show that UbiGate does not introduce any significant performance penalty in service discovery while enabling support of new services to be easily added.

The rest of this paper is organized as follows. Section 2 introduces our approach and the architecture of the UbiGate gateway. Section 3 assesses our prototype implementation of UbiGate, demonstrating its benefits in a pervasive environment. In Section 4, we present other approaches that were proposed as attempts to unify the mechanisms of service discovery. Finally, Section 5 concludes and presents future work. The most important concern in service discovery within pervasive computing is the interconnectivity among communication entities. Yet, the variety of protocols, because of their incompatibility, has made it difficult for one device to keep track of the availability of existing services. The underlying challenges are to enable service discovery using a unique protocol, and to improve reactivity, considering the volatility of pervasive environments. Our approach consists of building a gateway in charge of transforming service discovery messages into presence information for end-users. This gateway, called UbiGate, exchanges information with users through the SIP protocol, while managing various native discovery protocols, as illustrated in Figure 1.

OUR APPROACH

UbiGate aims at unifying service discovery protocols so as to make them transparent for users. Besides, it is pertinent to suppose that the appearance of new protocols will facilitate the advertisement of more services that may interest users. Therefore, the infrastructure has been designed to easily welcome new protocols. In UbiGate, protocols are plugged as modules packaging their capabilities.

Practically, UbiGate enables service discovery through SIP compliant devices. Using the SIP protocol, users can address discovery requests to the gateway which will manage all tasks for information gathering. SIP is thus the pillar of our approach. So as to better understand how it sustains the infratructure, an overview of the protocol is presented.

SIP background

SIP is originally a signaling protocol for Voice over IP (VoIP) and third generation mobile phones. It is standardized by the IETF and adopted by the ITU. 1 This protocol enables creating, modifying and terminating a communication between parties. Communications include audio/video communications, games, and instant messaging.

The SIP protocol is based on a client-server model. A SIP message can be sent or received by a client. A sent message is said to be outgoing; otherwise, it is said to be incoming.

SIP is a text-based protocol similar to other well-known protocols such as HTTP2 and RTSP. 3 A SIP message begins with a line indicating whether the message is a request (including a protocol method name) or a response (including a return code). A sequence of required and optional headers follows. Finally, a SIP message includes a body containing other information relevant to the message.

Logically, SIP is composed of three main entities: a registrar server to allow users to record their current location, a proxy server in charge of dispatching SIP messages, and a user agent required on each communication device to perform all SIP-related actions.

To support mobility, a user is assigned a SIP URI (Uniform Resource Identifier), which is a symbolic address, analogous to an e-mail address. When a SIP proxy receives a message for a local URI, it asks the local registrar server to translate the URI into contact information for a specific user agent. A user must thus inform the registrar server of the user agent at which he would like to receive messages.

Following the success of SIP, the IETF has produced many specifications related to presence and instant messaging with SIP. This set of specifications is known as SIMPLE 4 and covers topics ranging from protocols for subscription and publication to presence document formats.

UbiGate architecture

The UbiGate gateway aims at meeting the challenges of service discovery in pervasive computing by leveraging the SIP protocol. The architecture of UbiGate is based upon the three major components illustrated in Figure 2. The main issue that is addressed by the UbiGate gateway is the reliability of the communication channel between UbiGate and users. Indeed, service requesters rely on Ubi-Gate to perform discovery tasks and regularly convey back to them appropriate information. So as to ensure a reliable and efficient access to information, service attributes are stored in a centralized registry managed by the gateway. Indeed, because of the dynamicity of pervasive environments, distributed storage systems are hard to deploy and maintain efficiently. Considering the capabilities of SIP, UbiGate deploys an event framework in charge of managing a SIP presence server. This Presence Manager handles the storage of services as entries, and matches requesters' subscriptions with those entries.

In order to discover services, UbiGate uses a service broker. This broker handles the native discovery protocols supported by the gateway. The Communication Manager (i.e. the service broker), is in charge of activating service discovery sessions and recovering the attributes of services that are available in the environment. This component cooperates with the previous one for service information storage.

The Gateway Manager enforces the coordination between the components. It therefore enables them to conjointly perform discovery and notification tasks. This component also provides a control interface for the administration of UbiGate, allowing users to easily add new features.

Furthermore, since our approach aims at providing an infrastructure that is as much extensible as possible, Ubi-Gate's suported discovery protocols are integrated through modules. For a given discovery protocol, the module packaging its capacities can be plugged into the gateway and removed as well without impacting the behaviour of existing protocols. This design feature allows UbiGate to be easily integrated in the environment so as to limit the impact of changes for end-users. Typically, the gateway renders seamless the upgrades on discovery protocols as well as disruptions in environment settings. Indeed, service providers are not required to switch to or add new protocols in their service advertisement kernel. In point of fact, it is up to the gateway to adapt itself to the environment by integrating in its core new modules for supporting available discovery protocols.

Service Discovery using UbiGate

Service discovery through UbiGate can be split in two distinct activities: the discovery process activity which collects service information, and the notification activity during which information is delivered to service requesters.

UbiGate uses a Pull model for its discovery processes. Practically, all the native service discovery protocols managed by the gateway are regularly called upon to scan the environment so as to provide up-to-date information on available services. However, as far as service requesters are concerned, the discovery system is perceived as running a Push model that works in real time. Indeed, an advertisement is performed each time the context changes. This scheme allows UbiGate to unify service discovery mechanisms while improving the Push model. It is this notion of realtime notification that has encouraged the choice of SIP.

A user wishing to benefit from the capacities of UbiGate is required to register itself as a service requestor by subscribing to the presence of a service. The subscription sent to the gateway is performed by sending a SIP SUBSCRIBE request, including an URI5 that represents the targeted service. For the purpose of guaranteeing a uniform interaction with the gateway, we have introduced an ad-hoc format for encoding relevant subscription information into URIs. This format indicates the type of service, or/and the provider's name, and/or the preferences on the protocol through which the service is advertised. The service request format thus provided allows requesters to target a specific type of service on a given device while filtering the advertisement protocol. The latter option can be useful if the requester is sending requests on behalf of other entities that have strict constraints. As illustrated in the example of Figure 3, keywords such as Printer may be used so that communicating entities may understand the significance of the terms utilized.

To discover the full range of available services in its vicinity, the client can send beforehand a SUBSCRIBE [START_REF] Roach | Session Initiation Protocol (SIP)-Specific Event Notification[END_REF] request using the joker UbiAny. The gateway will then inform the client of all service types available through a list of keywords. Using them, the user can then properly transmit its subscription to the gateway. Service attributes requested by a user are described in a standard PIDF 6 XML document attached to SIP NOTIFY [START_REF] Roach | Session Initiation Protocol (SIP)-Specific Event Notification[END_REF] messages. To enable standard, unmodified, SIP user agents such as Windows Live Messenger to correctly process such a PIDF document, we use the optional note tag recommended by the RFC [START_REF] Sugano | Presence information data format (PIDF)[END_REF] for additional information, as illustrated in Figure 4. Each of these service descriptions contains only the essential information that characterizes a service: the service type, the address of the entity offering the service, and the native protocol used by the provider to announce the service.

Thus, the compact description proposed in Figure 4 presents a OBEX Object Push service advertised by the Bluetooth Service Discovery Protocol supported by a Bluetooth-enabled device which MAC address is '0C:19:73:1F:1C'. This type of description is provided when requesters specify in their requests the service type and the protocol. The presence server then supplies the provider's address as requested information.

However, it is noteworthy to mention that the flexibility of SIP can be exploited so as to deliver better descriptions of services. Indeed, besides the essential attributes listed previously, a user may be interested in knowing the name of the service provider, its state along with the possible functionalities it yields. To deliver such an extended description, we have designed a new Service Discovery Information Format that is an extension of PIDF. In this new format, instead of using the PIDF note tag, UbiGate fills a service tag with all attributes values. The new format thus allows a more enhanced description of services as illustrated in the example of Figure 5. The protocol tag in the Service Discovery Information Format contains a value that refers to the protocol used by the service provider to advertise the service. Thus, the value 'RFID' has been placed in the example of Figure 5 because an RFID tag has been attached to the printer Zeus, and it is the detection of this tag that has provided all the information. This extended description format has been designed to include all informations provided by Bluetooth Service Description Protocol as it was the experimental protocol used in UbiGate. Yet, the extensibility of the PIDF format guarantees future adaptations of the format to include any additional attribute proposed by any service discovery protocol. Moreover, the format proposed is flexible since the parser does not require all services to be completely detailed. Besides, note tags can be used when the description proposes more information that need to be conveyed to the requester.

Advanced Features. With UbiGate, a user can subscribe to any available services, without any concern for diversity of entities that provide those services, nor for the protocols they use for advertisement. However, usually, users require advanced discovery features based on service providers or more often on protocols used for service advertisement. Typically, if we consider a user that needs to discover, with his PDA, the services advertised by his desktop computer, it is important to allow him to subscribe to the computer rather than to all service types that UbiGate can find in the entire environment. Likewise, a user wishing to track all RFID tags that penetrate his environment should receive information related only to the protocol used for RFID. Thus, Ubi-Gate provides an optimisation of discovery requests. This optimisation consists of enabling a selection for the options mentionned above by including fields in subscription URIs.

ASSESSMENT

To assess our approach, we have developed a prototype implementation of UbiGate. This prototype includes about 10,500 lines of Java code, and relies on the SIP presence server 7 provided by NIST (National Institute for Standards and Technology).

For the purpose of presenting the contributions of Ubi-Gate, we have considered two case studies involving the Bluetooth technology and RFID which is increasingly adopted for tracking purpose.

First, we consider a building manager who needs to be informed in real time when items of great value cross specific gates of his building. These items have been associated with RFID tags. As an important constraint, the manager should be able to receive information on his desktop or laptop computer and any SIP compliant device. To enable this scenario, we have integrated in UbiGate support for RFID by developing a module, as described in Section 2.2. This module consists of about 200 lines of Java code. For our experiments, we used the ASPX RFID Kit (Icode 2) RW-310 with C-100 Converter.

In our second case study, we consider that a user with a Bluetooth-enabled device will likely scan the environment several times to discover available Bluetooth services. Since, inquiries in Bluetooth respect the Pull model, service discovery is no longer seamless nor effortless because it involves perpetual manipulation from users and/or permanent inquiries from devices. Using UbiGate, the user only has to subscribe to the presence of Bluetooth entities and their services. To do so, we have developed a Bluetooth module consisting of about 1,100 lines of Java code and based on the AvetanaBluetooth 8 stack.

As described in Section 2.3, creating well-formatted requests requires that users know their syntax, which can be a complex task. To overcome these limitations, we have developed a SIP User Agent that is compliant with UbiGate's new format for service description. It provides a powerful graphical interface simplifying the subscription steps. We have tested the presence server of UbiGate using SIP compliant standard clients such as SJPhone 9 and SIP Communicator 10 .

Within the same perspective, we have implemented an interface for UbiGate's administrators. They can use the interface to start pulling with a specific service discovery protocol even if there is no subscription yet. Thus, when requesters address their requests to UbiGate, the latter has already information to deliver. The interface also allows shutting down the gateway or disabling unused protocols.

Table 1 lists the sizes of the implementations for UbiGate's different components. It also includes the sizes of the two We now present preliminary results of a performance evaluation of our UbiGate implementation. For our experiments, the gateway has been deployed on an Intel Pentium 4 computer with CPU frequency of 3.06 GHz and 1 MB of memory. We have also setup a native application using the Avetana stack directly. Due to specificities of Bluetooth technology, only first inquiry sessions are relevant. Indeed, the service discovery process is then complete and the devices do not use any information from previous inquiries to accelerate information exchanges. The comparison in Table 2 shows that there is no significant overhead using Ubi-Gate compared to a direct discovery process using the Avetana Bluetooth stack. Indeed, over the 31 seconds required for bluetooth service discovery completion, UbiGate takes less than 1 second to reformat informations data and convey them through the network for display on the client side.

Native

UbiGate Delay 1 st inquiry ≈ 31 seconds ≈ 31.7 seconds

Table 2: Comparison of Bluetooth inquiry delays

We have also deployed Ubigate on an OSGI T M platform in order to facilitate its deployment, migration and the management of dependancies. For this, we have encapsulated UbiGate's components as OSGi bundles. The OSGi bundles for the modules enable dynamic upgrades of UbiGate, while presenting roughly the same size of code as the non-OSGi version. Our OSGI based implementation of UbiGate has been deployed and run successfully on the knopflerfish 11 OSGi platform. The OSGi bundles as well as the complete Java code of the UbiGate project are available at the project web page 12 .

Finally, it is noteworthy to indicate that UbiGate infrastructure allows new service discovery protocols (such as Jini, 11 Knopflerfish: www.knopflerfish.org 12 UbiGate: http://uuu.enseirb.fr/∼bissyand/UbiGate/ UPnP, ...) to be easily added. Further, as outlined in Figure 6, Ubigate enables IP and non-IP applications to discover each other transparently.

RELATED WORK

Over the years many research groups and industries have focused on designing solutions to cope with the heterogeneity of service discovery protocols. Typically, providing SDP interoperability consists of enabling applications to switch their current SDPs on the fly according to their networked environment. This is made possible through the use of an intermediate representation of SDPs paradigms (i.e an intermediary protocol) in order to abstract incompatibilities among SDPs to exclusively consider their similarities [START_REF] Becker | Base: A microbroker-based middleware for pervasive computing[END_REF][START_REF] Bromberg | INDISS: Interoperable discovery system for networked services[END_REF][START_REF] Costa | A Reconfigurable Component-based Middleware for networked Embedded Systems[END_REF][START_REF] Grace | A reflective framework for discovery and interaction in heterogeneous mobile environments[END_REF][START_REF] Loureiro | A flexible middleware for service provision over heterogeneous pervasive networks[END_REF][START_REF] Raverdy | A multi-protocol approach to service discovery and access in pervasive environments[END_REF]. So far, two approaches have emerged depending on how applications are either bound or unbound to this intermediary protocol [START_REF] Bromberg | Interoperability of service discovery protocols: Transparent versus explicit approaches[END_REF]. In an explicit approach, applications need to be explicitly designed to use a specific discovery API that translates, if required, the intermediary protocol to the SDP currently used according to the networked environment [START_REF] Raverdy | A multi-protocol approach to service discovery and access in pervasive environments[END_REF]. In a transparent approach applications are unaware of the translation process [START_REF] Bromberg | INDISS: Interoperable discovery system for networked services[END_REF]. While the latter offers seamless interoperability to legacy applications, the former enables the extending existing SDPs with advanced features. Our approach combines the strengths of the current approaches to provide a new SDP paradigm based on the concept of presence in order to provide a realtime discovery as required by pervasive environments. Following the example of the Amigo service architecture [START_REF] Georgantas | The amigo service architecture for the open networked home environment[END_REF], UbiGate allows the integration of heterogeneous technologies by establishing interoperability of SDPs through different mechanisms involving the deployment of an intermediate node.

CONCLUSION

The multiplicity of discovery protocols, the capacities of pervasive computers, and the expectations of users, have contributed to the rise of new challenges for service discovery. To meet these challenges, we have proposed a new approach for service discovery involving the use of a gateway and relying on a widespread protocol, SIP.

Following our approach, we have designed the UbiGate gateway in charge of transforming discovery information into presence information. UbiGate can then process subscriptions from SIP compliant devices and manage discovery tasks using adequate native service discovery protocols. Centered around the notion of presence, UbiGate enables a realtime availability of service information while hiding the heterogeneity of underlying protocols.

The suitability of the approach for pervasive computing has been illustrated through two case studies involving the use of the Bluetooth service discovery protocol and a protocol for RFID. We have then developed dedicated modules for UbiGate, thus demonstrating the ease for adding support of new protocols.

After tackling the issues in service discovery in this paper, future work will involve service delivery using the SIP protocol. Taking advantage of the peer-to-peer model of SIP and the portability aspect of UbiGate, service delivery can be dealt with efficiently.

Figure 1 :

 1 Figure 1: Interaction among pervasive computers, users and UbiGate

Figure 2 :

 2 Figure 2: UbiGate architecture

  sip:UbiAny.Printer.UbiAny@192.168.200.18

Figure 3 :

 3 Figure 3: Example of subscription URI

  <note> Bluetooth.OBEXObjectPush.0C:19:73:1F:1C </note>

Figure 4 :

 4 Figure 4: Compact description of a service

Figure 5 :

 5 Figure 5: Extract of the extended description of a service

Figure 6 :

 6 Figure 6: UbiGate service discovery stack (Example with Bluetooth SDP)

Table 1 :

 1 Package size of UbiGate components and modules

		Component/Module	Package size
		Gateway manager	4.8 Kb
	UbiGate	Communication manager	25.1 Kb
		SIP Presence Server	4.1 Mb
		Total size of UbiGate	4.13 Mb
	Modules	Bluetooth RFID	22.2 Kb 6.1 Kb

ITU: International Telecommunications Union.

HTTP: HyperText Transfer Protocol.

RTSP: Real Time Streaming Protocol.

SIMPLE: SIP for Instant Messaging and Presence Leveraging Extensions.

URI: Uniform Resource Identifier.

PIDF: Presence Information Data Format.

Acknowledgements

The authors would like to extend their gratitude to EN-SEIRB graduate students Nabila Ayadi, Jean Collas and Aurélie Goubault de Brugière who contributed major parts of the implementation of UbiGate's prototype. We further thank Dr. Julia Lawall and the anonymous reviewers for their useful comments.


Marie-Ange Raulet

Benjamin Ducharne

New mathematical

come

New mathematical approach using fractional derivatives to take into account excess losses in magnetic materials

Marie-Ange Raulet a , Benjamin Ducharne b , Daniel Guyomar Accurate prediction of magnetic circuit dynamical behaviour constitutes a fundamental problem for the design and the optimization of electromagnetic devices. The representation of dynamical effects due to eddy currents with a macroscopic model [1] is usually well admitted and adopted by researchers and electrical engineers. On the other hand, macroscopic modelling of microscopic dynamical effects due to wall motions which induce the so called excess losses, always exhibits interests of various research workers [2]- [6].

This paper proposes a new mathematical approach for the modelling of the dynamical microscopic effects. This approach is based on a fractional derivative formulation. This fractional term is used to describe the excess field contribution (1):

! ! " dt t B d t Hexc ) ( . ) ( # . ( 1
)
The fractional derivative model has already been used with success to model hysteresis on a large frequency bandwidth of ferroelectric ceramics (1mHz < ƒ < 10Hz) [7]. In the case of magnetic materials the fractional formulation is required to describe excess losses.

The new model is implemented and tested by considering a magnetic sample magnetized by a sine flux. An analytical study is carried out and gives results which are in a good agreement with the well known Bertotti's results related to the frequency relation of the excess energy [2].

Finally, comparisons between simulations and large number of experimental results performed on a SiFe sample magnetized by a varying frequency sine flux will be considered as validation of the new model.


Thibaut Chailloux

Laurent Morel

Fabien Sixdenier

Olivier Garrigues

Electromagnetic actuator to reduce vibration sources

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Electromagnetic actuator to reduce vibration sources

Thibaut Chailloux * , L. Morel * , F. Sixdenier * , O. Garrigues * * Université de Lyon, F69622, France ; Université Lyon 1, F69622, France ; CNRS UMR5005 AMPERE, 43, Bd du 11 Novembre 1918, Villeurbanne, F-69622, France E-mail : thibaut.chailloux@etu.univ-lyon1.fr , laurent.morel@univ-lyon1.fr , fabien.sixdenier@univ-lyon1.fr , garrigues@lagep.univ-lyon1.fr

In order to improve passenger comfort, a reduction of vibration sources in vehicles is being considered by the manufacturers. The vibrations can be compensated mainly by 4 ways : mechanical damping (passive system), hydraulic or pneumatic actuators, electromagnetic actuators [1], piezoelectric actuators.

As part of a study with Airbus, these solutions were compared: there is no universal solution, each should be used to compensate specific frequencies, amplitudes and forces, for a closed loop regulation.

In particular, Airbus wanted to compensate a vibration on an aircraft engine for which an electromagnetic actuator was recommended. (The specifications are : maximal force 6kN, frequency approximately 100Hz, movement amplitude 3µm, volume 10cm (length) x 10cm (width) x 20cm (height), temperature between -50 and +150°C).

Different topologies of electromagnetic actuators were compared and the chosen one was an electromagnet whose geometry has been optimized.

In regard of the 3D shape of the magnetic circuit, the choice of soft magnetic composite was quite natural because of its electromagnetic field isotropy.

In order to simulate the global system and estimate losses, we need an accurate model. [2], able to modelize dynamic hysteresis.

The aim of this study consists in finding a global behaviour model of the actuator that can give the temporal evolution of physical quantities and the different losses.

A first prototype that meets the main specifications has been built in the laboratory. The electromagnetic force, so as magnetic and electrical signals (voltage U, current I, magnetic flux !, magnetizing field H, magnetic induction B) generated by the actuator, and finally the actual losses, confirms that the soft magnetic composite is a good material candidate for this kind of actuators.

[ 1 ]

 1 Zupan et al./Actuator Classification and selection-The Development, Advanced Engineering Materials, Vol. 4, N°12 2002 [2] Dynamical Models for Eddy Current in Ferromagnetic Cores Introduced in an FE-Tuned Magnetic Equivalent Circuit of an Electromagnetic Relay, Sixdenier F., Raulet M.-A., Marion R., Goyet R., Clerc G., Allab F., IEEE Transactions on Magnetics 44, 6 (2008) 866-869

Figure 1 :

 1 Figure 1: General aspect of the actuator