Michael Emmi

Akash Lal
email: akashl@microsoft.com

Finding Non-Terminating Executions in Distributed Asynchronous Programs

Programming distributed and reactive asynchronous systems is complex due to the lack of synchronization between concurrently executing tasks, and arbitrary delay of message-based communication. As even simple programming mistakes have the capability to introduce divergent behavior, a key liveness property is eventual quiescence: for any finite number of external stimuli (e.g., client-generated events), only a finite number of internal messages are ever created.

In this work we propose a practical three-step reduction-based approach for detecting divergent executions in asynchronous programs. As a first step, we give a code-to-code translation reducing divergence of an asynchronous program P to completed state-reachability-i.e., reachability to a given state with no pending asynchronous tasks-of a polynomiallysized asynchronous program P . In the second step, we give a code-to-code translation under-approximating completed state-reachability of P by state-reachability of a polynomially-sized recursive sequential program P (K), for the given analysis parameter K ∈ N. Following Emmi et al.

[8]'s delay-bounding approach, P (K) encodes a subset of P 's, and thus of P 's, behaviors by limiting scheduling nondeterminism. As K is increased, more possibly divergent behaviors of P are considered, and in the limit as K approaches infinity, our reduction is complete for programs with finite data domains. As the final step we give the resulting state-reachability query to an off-the-shelf SMT-based sequential program verification tool. We demonstrate the feasibility of our approach by implementing a prototype analysis tool called Alive, which detects divergent executions in several hand-coded variations of textbook distributed algorithms. As far as we are aware, our easy-to-implement prototype is the first tool which automatically detects divergence for distributed and reactive asynchronous programs.

Introduction

The ever-increasing popularity of online commercial and social networks, along with proliferating mobile computing devices, promises to make distributed software an even more pervasive component of technological infrastructure. In a distributed program a network of physically separated asynchronous processors coordinate by sending and asynchronously receiving messages. Such systems are challenging to implement because of several uncertainties, including processor timings, message delays, and processor failures. Although simplifying mechanisms such as synchronizers and shared-memory simulation do exist [START_REF] Lynch | Distributed Algorithms[END_REF], they add significant runtime overhead which can be unacceptable in many situations.

Because of the inherit complexity in distributed asynchronous programming, even subtle design and programming mistakes have the capability to introduce erroneous or divergent behaviors, against which the usual reliability measures are much less effective. The great amount of nondeterminism in processor timings and message delays tends to make errors elusive and hard to reproduce in simulation and testing. The combinatorial explosion incurred by the vast number of processor interleavings and message-buffer contents tends to make formal verification techniques intractable. Though many distributed algorithms are proposed along with manual correctness proofs, key properties such as eventual quiescence-i.e., for any number of external stimuli such as client-generated events, only a finite number of internal network messages are ever createdremain difficult to ensure with automatic techniques. Practically speaking, such properties ensure the eventual construction of network spanning trees [START_REF] Lynch | Distributed Algorithms[END_REF], the eventual election of network leaders [START_REF] Svensson | A new leader election implementation[END_REF], and the eventual acceptance of network peer proposals, e.g., according to the Paxos protocol [START_REF] Lamport | The part-time parliament[END_REF].

In this work we develop an automatic technique to detect violations to eventual quiescence, i.e., executions of distributed systems for which a finite number of external stimuli result in an infinite number of internal messages. Our reduction-based approach works in three steps. First, we reduce the problem of finding nonterminating executions of a given (distributed) asynchronous program P to the problem of computing reachability in a polynomially-sized (distributed) asynchronous program P . This reduction is complete for programs with finite data domains, in the sense that an answer to the reachability query on P is a precise answer to the nontermination query on P . In the second step, we reduce reachability in P to reachability in a polynomially-sized recursive sequential program P -without explicitly encoding the concurrent behavior of P as data in P . This step is parameterized by an integer K ∈ N; for small K, P encodes few concurrent schedules of P ; as K is increased, P encodes and increasing number of concurrent reorderings, and in the limit as K approaches infinity, P codes all possible behaviors of P -and thus P . Finally, using existing sequential program verification tools, we check reachability in P : a positive result indicates a nonterminating execution in P , though the lack of nonterminating executions in P can only be concluded in the limit as K approaches infinity. Our technique supports fairness, in that we may consider only infinite executions in which no message is ignored forever.

We demonstrate the feasibility of our reduction-based approach by implementing a prototype analysis tool called Alive, which detects violations to eventual quiescence in several hand-coded variations to textbook distributed algorithms [START_REF] Lynch | Distributed Algorithms[END_REF]. Our relatively easy-to-implement prototype leverages existing SMT-based program verification tools [START_REF] Lal | Corral: A solver for reachability modulo theories[END_REF], and as far as we are aware, is the first tool which can automatically detect divergence in distributed asynchronous programs.

To begin in Section 2, we introduce a program model of distributed computation. In Section 3 we describe our reduction to sequential program analysis, and provide code-to-code translations which succinctly encode the reduction. Following in Section 4 we describe our experimental results in analyzing textbook distributed algorithms, and we conclude by discussing related work in Section 5.

Distributed Asynchronous Programs

We consider a distributed message-passing program model in which each processor is equipped with a procedure stack and an unordered buffer of pending messages. Initially all processors are idle. When an idle processor's message buffer is nonempty, some message is removed, and a message-dependent task is executed to completion. Each task executes essentially as a recursive sequential program, which besides accessing its own processor's global storage, can post messages to the buffers of any processor, including its own. When a task does complete, its processor again becomes idle, chooses a next pending message to remove, and so on. The distinction between messages and handling tasks is purely aesthetic, and we unify the two by supposing each message is a procedure-and-argument pair. Though in principle many message-passing systems, e.g., in Erlang and Scala, allow reading additional messages at any program point, we have observed that common practice is to read messages only upon completing a prior task [START_REF] Trottier-Hebert | Learn you some Erlang for great good[END_REF].

Our choice to model message-passing programs with unordered buffers has two important consequences. First, although some programming models do not ensure messages are received in the order they are sent, others do; our unordered buffer model should be seen as an abstraction of a model with faithful message queues, since ignoring message order allows behaviors infeasible in the queue-ordered model. Second, when message order is ignored, distributed executions are task-serializable-i.e., equivalent to executions where the tasks across all processors execute serially, one after the other. Intuitively this is true because (a) tasks of different processors access disjoint memory, and (b) message posting operations commute with each other. (Message posting operations do not commute when buffers are ordered.) To simulate a distributed system with a single processor we combine each processor's global storage, and ensure each processor's tasks access only their processor-indexed storage. Since serializability implies that single processor systems precisely simulate the behavior of distributed systems, we limit our discussion, without loss of generality, to single-processor asynchronous programs [START_REF] Sen | Model checking multithreaded programs with asynchronous atomic methods[END_REF]. Appendix A illustrates this mapping to a singleprocessor asynchronous program in more detail.

Program Syntax

Let Procs be a set of procedure names, Vals a set of values, Exprs a set of expressions, Pids a set of processor identifiers, and let T be a type. Figure 1 gives the grammar of asynchronous message-passing programs. We intentionally leave the syntax of expressions e unspecified, though we do insist Vals contains true and false, and Exprs contains Vals and the (nullary) choice operator . We say a program is finite-data when Vals is finite. Each program P declares a single global variable g and a procedure sequence, each p ∈ Procs having a single parameter l and top-level statement denoted s p ; as statements are built inductively by composition with control-flow statements, s p describes the entire body of p. The set of program statements s is denoted Stmts. Intuitively, a post p e statement is an asynchronous call to a procedure p with argument e. The assume e statement proceeds only when e evaluates to true, and this statement plays a role in disqualifying executions in our subsequent reductions of Section 3. The programming language we consider is simple, yet very expressive, since the syntax of types and expressions is left free, and we lose no generality by considering only single global and local variables. Appendix B lists several syntactic extensions which we use in the source-to-source translations of the subsequent sections, and which easily reduce to the syntax of our grammar.

Program Semantics

A (procedure) frame f = , s is a current valuation ∈ Vals to the procedurelocal variable l, along with a statement s ∈ Stmts to be executed. (Here s describes the entire body of a procedure p that remains to be executed, and is initially set to p's top-level statement s p ; we refer to initial procedure frames t = , s p as tasks, to distinguish the frames that populate task buffers.) The set of all frames is denoted Frames. A configuration c = g, w, m is a current valuation g ∈ Vals to the processor-global variable g, along with a procedure-frame stack w ∈ Frames * and a multiset m ∈ M[Frames] representing the pending-tasks buffer. The configuration c is called idle when w = ε, and completed when w = ε and m = ∅. The set of configurations is denoted Configs. The Post rule creates a new frame to execute the given procedure, and places the new frame in the pending-tasks buffer. The Complete rule returns from the final frame of a task, rendering the processor idle, and the Dispatch rule schedules a pending task on the idle processor.

An execution of a program P (from c 0 ) is a configuration sequence ξ = c 0 c 1 . . . such that c i → c i+1 for i ≥ 0; we say each configuration c i is reachable from c 0 . An initial condition ι = g 0 , 0 , p 0 is a global-variable valuation g 0 ∈ Vals, along with a local-variable valuation 0 ∈ Vals, and a procedure p 0 ∈ Procs. A configuration c = g 0 , 0 , s p0 , ∅ of a program P is called g 0 , 0 , p 0 -initial. An execution ξ = c 0 c 1 . . . is called infinitely-often idle when there exists an infinite set I ⊆ N such that for each i ∈ I, c i is idle.

Definition 1 (state-reachability). The (completed) state-reachability problem is to determine for an initial condition ι, global valuation g, and program P , whether there exists a (completed) g-valued configuration reachable in P from ι.

In this work we are interested in detecting non-terminating executions due to asynchrony, rather than the orthogonal problem of detecting whether each individual task may alone terminate. Our notion of non-termination thus considers only executions which return to idle configurations infinitely-often.

Definition 2 (non-termination). The non-termination problem is to determine for an initial condition ι and a program P , whether there exists an infinitelyoften idle execution of P from some ι-initial configuration.

Detecting Non-Termination

Though precise algorithms for detecting (fair) non-termination in finite-data asynchronous programs are known (see Ganty and Majumdar [START_REF] Ganty | Algorithmic verification of asynchronous programs[END_REF]), the fair nontermination problem is polynomial-time equivalent to reachability in Petri nets, which is an EXPSPACE-hard problem for which only non-primitive recursive algorithms are known. Though worst-case complexity is not necessarily an indication of feasibility on practically-occurring instances, here we are interested in leveraging existing tools designed for more tractable problems whose solutions can be used to incrementally under-approximate non-termination detection; i.e., where for a given analysis parameter k ∈ N we can efficiently detect non-termination from an interesting subset B k of program behaviors.

Our strategy is to reduce the problem of detecting non-terminating executions in asynchronous programs to that of completed state-reachability in asynchronous programs. We perform this step using the code-to-code translation of Section 3.1, and in Section 3.2 we consider extensions to handle fairness. Then, in the second step of Section 3.3, we apply an incrementally underapproximating reduction from state-reachability in asynchronous programs to state-reachability in sequential program [START_REF] Emmi | Delay-bounded scheduling[END_REF][START_REF] Bouajjani | On sequentializing concurrent programs[END_REF], and discharge the resulting program analysis problem using existing sequential analysis tools.

Reduction from Non-Termination to Reachability

In the first step of our reduction, we use the fact that every infinite execution eventually passes through two configurations c 1 , and then c 2 , such that every possible execution from c 1 is also possible from c 2 ; e.g., when c 1 and c 2 are idle configurations with the same global valuation in which all tasks pending at c 1 are also pending at c 2 . Formally, given two configurations c 1 = g 1 , w 1 , m 1 and c 2 = g 2 , w 2 , m 2 we define the order c 1 c 2 to hold when g 1 = g 2 , w 1 = w 2 , and m 1 ⊆ m 2 . 3 An execution c 0 c 1 . . . is called periodic when c i c j for two idle configurations c i and c j such that i < j. 4 The following lemma essentially exploits the fact that is a well-quasi-ordering on idle configurations. Lemma 1. A finite-data program P has an infinitely-often idle execution from ι if and only if P has a periodic execution from ι.

Proof. Suppose c 0 c 1 . . . is the sequence of idle configurations in an infinitely-often idle execution ξ. As the subset order ⊆ on multisets is a well-quasi order, and the domain Vals of global variables is finite, is a well-quasi order on idle configurations. Thus there exists i < j such that c i c j , so ξ is also periodic.

Supposing ξ = c 0 c 1 . . . is a periodic execution from ι, there exists idle configurations c i and c j of ξ such that i < j and c i c j ; let c i = g i , ε, m i and c j = g j , ε, m j . Since g i = g j and m i ⊆ m j , by definition of , the sequence of execution steps between c i and c j is also enabled from configuration c j -we may simply ignore the extra tasks m j \ m i pending in c j . For any k, l ∈ N and task buffer m ∈ M[Frames] such that k < l < |ξ|, let ξ m k,l be the sequence of configurations c k c k+1 . . . c l-1 of ξ, each with additional pending tasks m. Furthermore, let k • m be the multiset union of k copies of m. Letting m = m j \ m i , then ξ 0,i ξ i,j ξ m i,j ξ 2m i,j ξ 3m i,j . . . is an infinitely-often idle execution from ι which periodically repeats the same transitions used to construct ξ between c i and c j .

We reduce the detection of periodic executions to completed state reachability in asynchronous programs. Essentially, such a reduction must determine multiset inclusion between the unbounded task buffers at two idle configurations; i.e., for some idle configuration c i = g i , ε, m i reachable in an execution c 0 c 1 . . ., there exists j > i such that c j = g j , ε, m j with g i = g j and m i ⊆ m j . As the set m i of pending tasks at c i is unbounded, any reduction cannot hope to store arbitrary m i for later comparison with m j using finite-domain program variables.

Our reduction determines the correspondence between unbounded task buffers in the source program using only finite-domain program variables by leveraging the task buffers of the target program. For each instance of a task t which is pending in c i , we post an additional task pro(t) when t is posted; for each task t pending in c j , we either post an additional task anti(t) instead of t, or we post nothing, to handle the case where t is never dispatched. We then check that for each executed pro(t) a matching anti(t) is also executed, and that at some point no pro(t) nor anti(t) tasks are pending. By considering executions which alternate between tasks of {pro(t) : t ∈ m i } and {anti(t) : t ∈ m j }-where m j ⊆ m j such that m j \ m j correspond to the dropped tasks-we can ensure each instance of an m i task has a corresponding instance in m j , storing only the last encountered pro(t) task, for t ∈ m i .

Figure 3 lists our code-to-code translation ((P )) nt reducing non-termination in an asynchronous program P to completed state reachability in the asynchronous program ((P )) nt . Besides the auxiliary variable last used to store the last encountered pro(t) task, for t ∈ m i , we introduce Boolean variables repeated, to signal whether m i is non-empty, and turn, to signal whether an anti(t) task has been executed since the last executed pro(t) task. We also divide the execution of tasks into two phases by introducing a task-local Boolean variable period. The first phase (!period) corresponds to the execution c 0 c 1 . . . c i , while the second phase (period) corresponds to c i+1 c i+2 . . . c j . Initially pending tasks occur in the first non-period phase. Then each time a new task t is posted, a non-deterministic choice is made for whether t will execute in the non-period phase, in the period phase, or never.

Finally, to determine which finite asynchronous executions prove the existence of infinite asynchronous executions, we define the predicate ϕ nt over initial conditions ι and configuration c as

ϕ nt (ι, c) def =        true
when ¬repeated(ι) and turn(ι) and repeated(c) and turn(c) and

G[0](c) = G[1](ι) = G[1](c) false otherwise,
along with the mapping ϑ nt which projects the initial conditions of ((P )) nt to those of P , as ϑ nt ( g, , p ) def = g , , p when g(g ) = G[0](g), l( ) = l( ), and p = p. Essentially, in any completed configuration c reachable from ι satisfying ϕ nt (ι, c), we know that some task has executed during the period (since repeated evaluates to true), and that for each task pending at the beginning of the period, an identical task is pending at the end of the period (since turn evaluates to true, and there are no pending tasks in c). Finally, the conditions on the global variable G ensure that the beginning and end of each period reach the same global valuation.

Lemma 2. A finite-data program P has an infinitely-often idle execution from ι 0 if and only if a completed configuration c is reachable in ((P )) nt from some ι such that ϕ nt (ι, c) = true and ϑ nt (ι) = ι 0 .

Proof. For the forward direction, by Lemma 1, P also has a periodic execution ξ = ξ 0,i ξ i,j ξ j,ω from ι 0 -where ξ k,l def = c k c k+1 . . . c l-1 for k < l < |ξ|-and c i c j for idle configurations c i = g, ε, m 1 and c j = g, ε, m 2 . We build an execution ξ = ξ 0,i ξ i,j ξ match of ((P )) nt such that

-the configurations c k of ξ 0,i correspond to configurations c k of ξ 0,i , with g(c k ) = G[0](c k ), G[1](c k ) = g, -the configurations c k of ξ i,j correspond to configurations c k of ξ i,j , with g(c k ) = G[1](c k ) and G[0](c k ) = g, - the
pending tasks of each configuration c k of ξ 0,j , excluding pro and anti tasks, are contained within those of c k , the local valuations of each configuration c k of ξ 0,i (resp., of ξ i,j ) match those of c k , except period evaluates to 0 (resp., to 1) in every frame of c k , and the sequence ξ match alternately executes pro and anti tasks such that each pro(t) task is followed by a matching anti(t) task.

It follows that we can construct such a ξ which reaches a completed configuration c from some ι such that ϕ nt (ι, c), ϑ nt (ι) = ι 0 , and

G[0](c) = G[1](c) = g.
For the backward direction, the reachability of a completed configuration c of ((P )) nt from ι such that ϕ nt (ι, c) implies that there exists a periodic execution ξ = c 0 c 1 . . . of P ; in particular, there exist configurations c i c j of ξ with i < j, and which have the global valuations g(c i ) = g(c j ) = G[0](c) = G[1](c) reached at the end of each period of ((P )) nt 's execution, and the set of pending tasks m in c i are those second-period tasks posted by ((P )) nt from first-period tasks. Since the set of tasks posted and pending by the end of the second period must contain m-otherwise unexecutable pro tasks would remain pending-we can construct a run where the pending tasks of c j contain the pending tasks of c i , and so P has a periodic execution. By Lemma 1 we conclude that P also has an infinitely-often idle execution.

Ensuring Scheduling Fairness

In many classes of asynchronous systems there are (at least) two sensible notions of scheduling fairness against which to determine liveness properties: an infinite execution is called strongly-fair if every infinitely-often enabled transition is fired infinitely often, and weakly-fair if every infinitely-often continuously enabled transition is fired infinitely often. In our setting where asynchronous tasks execute serially from a task buffer, weak fairness becomes irrelevant; while one task executes no other transitions are enabled, and when idle (i.e., while no tasks are executing), all pending tasks become enabled. Furthermore once a task is posted, it becomes pending, and it is thus enabled in all subsequent idle configurations until dispatched. We thus define fairness according to what is normally referred to as strong fairness: an execution is fair when each infinitely-often posted task is infinitely-often dispatched.

To extend our reduction so that only fair infinite executions are considered we make two alterations to the translation of Figure 3. First, on Line 36 we replace skip with assume period; this ensures participation of all tasks pending at the beginning of each period. Second, we add auxiliary state to ensure at least one instance of each task posted during the period is dispatched. This can be encoded in various ways; for instance, we can add two arrays dropped and dispatched of index type Procs × Vals and element type B that indicate whether each task has been dropped/dispatched during the period phase (i.e., where the local variable period evaluates to true). Initially dropped[t] = dispatched[t] = false for all t ∈ Procs × Vals. Each time a post to task t is dropped during the period phase (i.e., Line 36) we set dropped[t] to true, and each time task t is executed during the period phase (i.e., Line 38 when period is true) we set dispatched[t] to true. (Note that we need not consider the non-post of t on Line 34 as dropped, since t is necessarily dispatched during the period phase; otherwise there would remain a pending anti(t) task.) Finally, we add to our reachability query the predicate ∀t.dropped[t]⇒dispatched[t], thus ensuring that when all asynchronous tasks have completed the only dropped tasks have been dispatched during the period.

Alternatively, we may also encode this fairness check by posting auxiliary dropped and dispatched tasks to the task buffer, in place of using the dropped and dispatched arrays. Essentially for each task t dropped during the period phase on Line 36 we add post dropped(t), and for each task t posted into the period phase we add post dispatched(t). Then, using a single additional variable of type Procs × Vals we ensure that for every executed dropped(t) task some dispatched(t) task also executes; a single variable suffices for this check because we may consider only schedules where all dropped(t) and dispatch(t) tasks execute contiguously for each t.

Delay-Bounded Reachability

Following the reduction from (fair) nontermination, we are faced with a highlycomplex problem: determining completed state-reachability in finite-data programs is polynomial-time equivalent to computing exact reachability in Petri nets (i.e., such that all places representing pending tasks are empty), or alternatively in vector addition systems (i.e., such that all vector components counting pending tasks are zero). Though these problems are known to be decidable, there is no known primitive-recursive upper complexity bound.

Rather than dealing with such difficult problems, our strategy is to consider only a restricted yet interesting set of actual program behaviors. Following Emmi et al. [START_REF] Emmi | Delay-bounded scheduling[END_REF]'s delay-bounding scheme, we equip some deterministic task scheduler with the ability to deviate from its deterministic schedule only a bounded number of times (per task). As this development is very similar to Emmi et al. [START_REF] Emmi | Delay-bounded scheduling[END_REF]'s, we 1 // translation of var g: T 2 var g: T 3 var G[K]: T refer the interested reader there. We recall in Figure 4 the essential delay-bounded asynchronous to sequential translation.

To determine which executions of the sequential program ((P ))

K db prove the existence of a valid asynchronous execution, we define the predicate ϕ db over initial conditions ι and configuration c as

ϕ db (ι, c) =      true when G[0](ι) = g(c) and ∀i ∈ N.0 < i < K ⇒ G[i](ι) = G[i -1](c) false otherwise,
along with the mapping ϑ db from initial conditions of ((P ))

K db to those of P as ϑ db ( g, , p ) def = g , , p when g(g ) = g(g), l( ) = l( ), and p = p. Essentially, in any completed configuration c reachable from ι satisfying ϕ db (ι, c), we know that the initially pending task returned with the shared global valuation G[0](ι) resumed by the first-round tasks, and that the last (i-1) round task, for 0 < i < K, returned with the shared global valuation G[i](c) resumed by the first i round task. The following lemma follows from Emmi et al. [START_REF] Emmi | Delay-bounded scheduling[END_REF]. Lemma 3. A valuation g is reachable in some completed configuration of a program P from ι 0 if some g-valued completed configuration c is reachable in ((P )) K db from some ι, such that ϕ db (ι, c) = true and ϑ db (ι) = ι 0 , for some K ∈ N.

Experience

We have implemented a prototype analysis tool called Alive. Our tool takes as input distributed asynchronous programs written in a variation of the Boogie language [START_REF] Barnett | Boogie: An intermediate verification language[END_REF] in which message posting is encoded with specially-annotated procedure calls. Given a possibly non-terminating input program P , Alive translates P into another asynchronous program P (according to the translation of Sections 3.1 and 3.2) that may violate a particular assertion if and only if P has a (fair) non-terminating execution. Then Alive passes P and a bounding parameter K ∈ N to our AsyncChecker delay-bounded asynchronous program analysis tool [START_REF] Emmi | Asynchronous programs with prioritized task-buffers[END_REF] which attempts to determine whether the assertion can be violated (in an execution using at most K delay operations, per task). AsyncChecker essentially performs a variation of our delay-bounded translation of Section 3.3which results in a sequential Boogie program-and hands the resulting program P to the Corral SMT-based bounded model checker [START_REF] Lal | Corral: A solver for reachability modulo theories[END_REF] to detect assertion violations.

Our implementation is able to find (fair) non-terminating executions in several toy examples, and handed-coded implementations of textbook distributed algorithms [START_REF] Lynch | Distributed Algorithms[END_REF]; the source code of our examples can be found online [START_REF] Emmi | Finding non-terminating executions in distributed asynchronous programs[END_REF]. Figure 5 summarizes our experiments on three families of examples which we discuss below: the PingPong family of toy examples, the SpanningTree family of textbook examples, and variations on Lamport's Paxos algorithm [START_REF] Lamport | The part-time parliament[END_REF]. For each family, Figure 5 lists both "buggy" variations (i.e., those with infinite executions) and "correct" variations (those without infinite executions-at least up to the given delay bound). In each case the delay bound is given by K, and a recursion bound is given by N ; our back-end bounded model checker Corral only explores executions in which the procedure stack never contains more than N frames of any procedure, for a given recursion bound N ∈ N. Note that our implementation is a simple unoptimized prototype; the running times are simply listed as a validation that our reduction is feasible.

PingPong

As a simple example of a non-terminating asynchronous program, consider the PingPong program of Figure 6. Initially the Main procedure initializes the Boolean variable x to false and posts asynchronous calls to Ping and Pong. When Ping executes and x is false, then Ping posts a subsequent call to Ping, and sets x to true; otherwise Ping simply returns. Similarly, when Pong executes and x is true, then Pong posts a subsequent call to Pong, and sets x to false; otherwise Pong simply returns. This program has exactly one non-terminating execution: that where the pending instances to Ping and Pong execute in alternation. This execution is periodic, as the configuration where x=false, and both Ping and Pong have a single pending instance, is encountered infinitely often.

Figure 6a depicts an execution of the program resulting from our translation (Section 3.1) of the PingPong program. Following our translation, the Main procedure takes the branch of Line 28 in Figure 3, posting both pro(Ping) and Ping, then both pro(Pong) and Pong. Without using any delay operations, the scheduler encoded by AsyncChecker executes the posted tasks in depthfirst order over the task-creation tree [START_REF] Emmi | Delay-bounded scheduling[END_REF][START_REF] Emmi | Asynchronous programs with prioritized task-buffers[END_REF]. Thus following Main, pro(Ping) executes, then Ping, followed by anti(Ping). Subsequently, pro(Pong), Pong, and anti(Pong) execute, in that order. Luckily this execution provides a witness to nontermination without spending a single delay.

Our experiments include several variations of this example. The -mod2 and -mod3 variations add an integer variable i which is incremented (modulo 2, resp., 3) by each call of Ping. The addition of this counter complicates the search for a repeated configuration, since besides the global variable x and pending tasks Ping and Pong, the value of i must also match in the repeating configuration. This addition also increases the number of delay operations required to discover an infinite execution, as the depth-first task scheduler without delaying considers only executions where all Ping tasks execute before all Pong tasks-see Figure 6b; since, for instance, modulo 2 incrementation requires two of each Ping and Pong tasks to return to a repeating configuration (i.e., with i=0), the second Ping task must delay in order to occur after the first Pong task. In the -1md and -2md variations, we reduce the budget of task delaying, and observe that indeed the additional delay budgets are required to witness nonterminating executions. The PingPongPung variation is an even more intricate variation in which each task (i.e., Ping, Pong, or Pung) posts a different task.

SpanningTree

In Figure 7 we consider two examples of distributed algorithms taken from the textbook of Lynch [START_REF] Lynch | Distributed Algorithms[END_REF], and modified to introduce nonterminating executions. Essentially, SpanningTree attempts to compute a spanning tree for an arbitrary network by building a parent relation from message broadcasts. When the parent link is established asynchronously there exist (unfair) executions in which nodes cyclically propagate their search messages without ever establishing the parent relation. The BellmanFord algorithm is a generalization of SpanningTree in which links between nodes have weights; the algorithm attempts to establish a spanning tree in which each node is connected by a minimal-weight path. Our injection of a bug demonstrates that even the most trivial of programming errors (e.g., typing ≤ rather than <) can introduce fair nonterminating executions.

Alive automatically discovers these nonterminating executions for an arbitrary, unspecified network.

Paxos

Lamport's Paxos algorithm [START_REF] Lamport | The part-time parliament[END_REF] provides a two-phase protocol for collaboratively choosing a (numeric) value from a set of values proposed by various nodes in a network; Figure 8 lists a basic variation of the algorithm. Initially a set of proposers choose a unique value to propose, and broadcast their intention to the set of acceptors via the prepare message. Each acceptor then decides whether to support the proposed value, depending on whether or not a higher proposal has already been seen. When a proposal_OK message is received, the proposer checks whether a majority has been achieved, and if so broadcasts an accept message.

If in the meantime the acceptors have not encountered a higher proposal, they agree on the given proposal by setting accepted on Line 46. Even in fair executions, divergent behavior can arise from several places. As in the program of Figure 8, the proposers may periodically post higher proposals in case their initial proposal is not answered within a timeout (Line 12), when NOTIFY_DECLINED is false. Then even an individual proposer may repeatedly propose new values just before receiving the acceptors' proposal_OK messages. The acceptors, in turn, may continue to increment their prepared values, such that previously agreed proposals will no longer be accepted (see the condition on Line 40). Even preventing such behavior by assuming the proposers only submit new proposals upon the reception of declined messages (i.e., suppose NOTIFY_DECLINED is true), fair nonterminating executions may still arise by competition between two or more proposers; for instance where two proposers continuously outbid the other before either's proposal has been accepted.

Since each subsequent proposal in the Paxos algorithm proposed an increasingly large number, strictly speaking our detection algorithm will not discover such nonterminating executions, since the same values of proposal and prepared will not be encountered twice. Essentially we must extend our well-quasi-ordering of Section 3.1 by relaxing the equality on global state valuations to a well-quasiordering which is compatible with the program's transition relation. For the purpose of our experiments, we have encoded manually such an order for our variations on the Paxos algorithm; the order relates global valuations g 1 g 2 when there exists some δ ∈ N such that the values of proposal for proposing processes, and prepared for accepting processes, in g 1 and g 2 uniformly increase by δ, and all other variables in g 1 and g 2 are equal. With this small manual effort, Alive is able to discover the "individual" nonterminating execution described above, and while Alive can also detect the "competing" nonterminating execution in theory, AsyncChecker times out on the reachability check after 30 minutes.

Related Work

Contrary to much work on sequential program (non)termination detection [START_REF] Cook | Termination proofs for systems code[END_REF][START_REF] Gupta | Proving non-termination[END_REF], less attention has been paid to concurrent programs, where nontermination can arise from asynchronous interaction rather than diverging data values. Though both Cook et al. [START_REF] Cook | Proving thread termination[END_REF] and Popeea and Rybalchenko [START_REF] Popeea | Compositional termination proofs for multithreaded programs[END_REF] have proposed techniques to prove termination in multithreaded programs, failure to prove termination does not generally indicate the existence of nonterminating executions. In very recent work, Atig et al. [START_REF] Atig | Detecting fair nontermination in multithreaded programs[END_REF] suggest compositional nontermination detection for multithreaded programs based on bounded context-switch; their technique detects infinite executions between a group of interfering, and each non-terminating, threads. Our approach is orthogonal, as we detect infinite executions in which every task terminates; nontermination arises from the never-ending creation of new tasks. Technically, while Atig et al. [START_REF] Atig | Detecting fair nontermination in multithreaded programs[END_REF] explore the behaviors between statically-known threads, our problem is to detect the repetition of an unbounded set of dynamically-created tasks.

Our reduction-based technique follows a recent trend of compositional translations to sequential program analysis by considering bounded program behaviors. Based on the notion of bounded context-switch [START_REF] Qadeer | Context-bounded model checking of concurrent software[END_REF], Lal and Reps [START_REF] Lal | Reducing concurrent analysis under a context bound to sequential analysis[END_REF] proposed

Conclusion

We have proposed a practical reduction-based algorithm for detecting divergent executions in distributed asynchronous programs. By incrementally increasing possible task reordering, our approach explores an increasing number of possiblydivergent behaviors with increasing analysis cost, and any possibly-divergent behavior is considered at some cost. By reducing divergence of distributed asynchronous programs to assertion violation in sequential programs, our approach leverages efficient off-the-shelf sequential program analysis tools. Using our prototype tool, Alive, we demonstrate that the approach is able to find divergent executions in modified versions of typical textbook distributed algorithms. The translated asynchronous program executes on a single processor and has a single unordered buffer of pending tasks. Its global memory is a union of the memory of each machine. The main procedure initializes all machines and starts the client. The sending of a message is simply a post to the task buffer of the asynchronous program. It is easy to see that the two programs have the same semantics.

B Syntactic Extensions Used in Our Code Translations

The following syntactic extensions are reducible to the original program syntax of Section 2.1. Here we freely assume the existence of various type-and expressionconstructors. This does not present a problem since our program semantics does not restrict the language of types nor expressions.

Multiple types. Multiple type labels T 1 , . . . , T j can be encoded by systematically replacing each T i with the sum-type T = j i=1 T i . This allows local and global variables with distinct types.

Multiple variables. Additional variables x 1 : T 1 , ..., x j : T j can be encoded with a single record-typed variable x: T , where T is the record type { f 1 : T 1 , ..., f j : T j } and all occurrences of x i are replaced by x.f i . When combined with the extension allowing multiple types, this allows each procedure to declare any number and type of local variable parameters, distinct from the number and type of global variables.

Local variable declarations. Additional (non-parameter) local variable declarations var l': T to a procedure p can be encoded by adding l' to the list of parameters, and systematically adding an initialization expression (e.g., the choice expression , or false) to the corresponding position in the list of arguments at each call site of p to ensure that l' begins correctly (un)initialized.

Unused values. Call assignments call x := p e, where x is not subsequently used, can be written as call _ := p e, where _: T is an additional unread local variable, or simpler yet as call p e.

Unused branches. if e then s else skip is abbreviated by if e then s. Let bindings. Let bindings of the form let x: T = e in can be encoded by declaring x as a local variable var x: T immediately followed by an assignment x := e. This construct is used to explicate that the value of x remains constant once initialized. The binding let x: T in is encoded by the binding let x: T = in where is the choice expression.

Arrays. Finite arrays with j elements of type T can be encoded as records of type { f 1 : T , ..., f j : T }, where f 1 . . . f j are fresh names. Occurrences of terms a[i] are replaced by a.f i , and array-expressions [ e 1 , ..., e j ] are replaced by record-expressions { f 1 = e 1 , ..., f j = e j }. procedure-frame stack. The Return rule removes the top-most procedure frame from the stack, and substitutes the valuation of the return expression e into the assignment x := left below by the matching call statement. Note that the transition relation → S is non-deterministic, since the evaluation of an expression e can result in an arbitrary set of possible values.

C Sequential Program Semantics

PFig. 1 .Fig. 2 .

 12 Fig.1. The grammar of asynchronous message-passing programs P . Here T is an unspecified type, and e and p range over expressions and procedure names.

Figure 2

 2 Figure 2 defines the transition relation → for the asynchronous behavior. (The transitions for the sequential statements are standard, and listed in Appendix C.)The Post rule creates a new frame to execute the given procedure, and places the new frame in the pending-tasks buffer. The Complete rule returns from the final frame of a task, rendering the processor idle, and the Dispatch rule schedules a pending task on the idle processor.An execution of a program P (from c 0 ) is a configuration sequence ξ = c 0 c 1 . . . such that c i → c i+1 for i ≥ 0; we say each configuration c i is reachable from c 0 . An initial condition ι = g 0 , 0 , p 0 is a global-variable valuation g 0 ∈ Vals, along with a local-variable valuation 0 ∈ Vals, and a procedure p 0 ∈ Procs. A configuration c = g 0 , 0 , s p0 , ∅ of a program P is called g 0 , 0 , p 0 -initial. An execution ξ = c 0 c 1 . . . is called infinitely-often idle when there exists an infinite set I ⊆ N such that for each i ∈ I, c i is idle.

Fig. 3 .

 3 Fig. 3. The translation ((P )) nt of an asynchronous program P .

4 5 / 8 9Fig. 4 .

 584 Fig. 4. The K delay sequential translation ((P )) K db of an asynchronous program P .

Fig. 6 .

 6 Fig. 6. The PingPong program, along with asynchronous executions of the translations ((PingPong)) nt (a) and ((PingPong-mod2)) nt (b). Task order is indicated by numeric prefixes; the dotted line indicates delaying.

Fig. 9 .

 9 Fig. 9. (a) A generic distributed program with N machines. (b) A semantically equivalent shared-memory asynchronous program with a single task buffer.

Increment.

  Increment operations x++ are encoded as x := x + 1.

For expressions without program

  variables, we assume the existence of an evaluation function • e : Exprs → ℘(Vals) such that e = Vals. For convenience we define e(g, ) def = e[g/g, /l] e to evaluate the expression e in a global valuation g by substituting the current values for variables g and l. As these are the only program variables, the substituted expression e[g/g, . . .] has no free variables. To further reduce clutter in the operational program semantics, we introduce a notion of context. A statement context S is a term derived from the grammar S ::= | S; s, where s ∈ Stmts. We write S[s] for the statement obtained by substituting a statement s for the unique occurrence of in S. Intuitively, a context filled with s, e.g., S[s], indicates that s is the next statement to execute in the statement sequence S[s]. Similarly, a configuration context C = g, , S w, m is a configuration whose top-most frame's statement is replaced with a statement context, and we write C[s] to denote the configuration g, , S[s] w, m . When e is an expression, we abbreviate e(C[skip]) by e(C).

Figure 10 definesFig. 10 .

 1010 Fig. 10. The transition relation → for the standard sequential program statements.

Here ⊆ is the multiset subset relation.

As our definition of only relates configurations with equal global valuations, our notion of periodic is only complete for finite-data programs.

* Supported by

a reduction from detecting safety violations in multithreaded programs (with a finite number of statically-known threads) to detecting safety violations in sequential programs; shortly after La Torre et al. [START_REF] La Torre | Model-checking parameterized concurrent programs using linear interfaces[END_REF] extended this result to handle an arbitrary number of parametric threads, which was further extended by Emmi et al. [START_REF] Emmi | Delay-bounded scheduling[END_REF] to handle dynamic thread creation-including the case of task-buffer based "asynchronous programs" [START_REF] Sen | Model checking multithreaded programs with asynchronous atomic methods[END_REF]. More recently Bouajjani and Emmi [START_REF] Bouajjani | Bounded phase analysis of message-passing programs[END_REF] proposed a reduction from safety violations in distributed asynchronous programs with ordered message queues. Thus far, only the recent (yet orthogonal-see above) work of Atig et al. [START_REF] Atig | Detecting fair nontermination in multithreaded programs[END_REF] considers liveness properties such as nontermination.

Finally, although reductions from fair nontermination of task-buffer based finite-data asynchronous programs (alternatively, Petri nets) are known-e.g., by encoding into Petri net path logic formalae [START_REF] Ganty | Algorithmic verification of asynchronous programs[END_REF]-our encoding into asynchronous programs is original, and takes advantage of existing program analysis tools with efficient under-approximating exploration strategies. Technically, Ganty and Majumdar [START_REF] Ganty | Algorithmic verification of asynchronous programs[END_REF]'s encoding uses constraints on marking-valued variables to ensure that each task pending at the beginning of a repeating period is re-posted and pending at the period's end; a path-logic solver must then determine satisfiability under those constraints. Our encoding handles the matching of pre-and postperiod pending tasks directly; we pose an asynchronous program reachability query on a program whose additional tasks block executions in which pre-and post-period tasks cannot be matched.


Farah Moety

Marc Ibrahim

Samer Lahoud

Kinda Khawam

Distributed Heuristic Algorithms for RAT Selection in Wireless Heterogeneous Networks

In wireless heterogeneous networks, one of the most challenging problems is Radio Access Technology (RAT) selection that must be designed to avoid resource wastage. In this paper we adopt a hybrid model for RAT selection where the system allocates the downlink traffic between two different technologies in order to enhance global performance. We study the case of an integrated hybrid Wireless Local Area Network environment where the challenge we face is the high computational complexity necessary to obtain the global optimal solution. Therefore, we propose four distributed heuristic algorithms for RAT selection, where two of them are based on the distance between the user and the access points (APs), namely, distance based and probabilistic distance based algorithms. While the two others schemes are based on the peak rate that each user receives from these APs (peak rate based and probabilistic peak rate based algorithms). Results show that the proposed algorithms give efficient results compared to the optimal one depending on the spatial users distribution. Moreover these algorithms have a low computational complexity which makes them more advantageous compared to the optimal scheme in presence of a large number of users.

I. INTRODUCTION

The next wireless generation is characterized by the coexistence of a variety of wireless access technologies such as Universal Mobile Telecommunications System (UMTS), Wireless Local Area Network (WLAN), Worldwide Interoperability for Microwave Access (WiMAX), Long Term Evolution (LTE), etc. These different RATs can be integrated together to complement each other in terms of coverage, bandwidth, services, and mobility support, thus providing benefits for both operator and users. The user connected to the "best network" benefits from a good quality of service (QoS) at low cost. The operator in turn can better exploit its radio resources, satisfy more customers and offer new services. In this heterogeneous context, a new challenging problem of RAT selection arises. RAT selection is the functionality devoted to decide to which RAT a given user will be associated, where the challenge is to optimize for example the QoS, the throughput, or the cost of the network. In [START_REF] Perez-Romero | Resource Management for Heterogeneous Wireless Access Networks, Heterogeneous Wireless Access Networks Architectures and Protocols[END_REF], three categories of RAT selection schemes are considered based on the parameters taken into account in decision making, namely, service-based, load balancing-based, and interferencebased algorithms. A service-based RAT selection scheme uses a direct mapping between services and RATs [START_REF] Perez-Romero | Radio Resource Management Strategies in UMTS[END_REF]. A load balancing strategy distributes the load among different RATs as evenly as possible [START_REF] Tolli | Performance Evaluation of Common Radio Resource Management (CRRM)[END_REF]; specifically, the selected RAT will be the one having the lowest load. An interference-based strategy takes transmitted power level, propagation conditions, interference, and specific network characteristics into consideration [START_REF] Perez-Romero | Policy-based Initial RAT Selection algorithms in Heterogeneous Networks[END_REF]. Different approaches for RAT selection were studied such as individual, global and hybrid approaches. In the individual approach, each user selfishly strives to improve its own performance; authors in [START_REF] Kumar | User-Network Association in an 802.11 WLAN & 3G UMTS Hybrid Cell: Individual Optimality[END_REF] studied this approach in an integrated WLAN and UMTS hybrid cell where users decide to join one of the two RATs so that their average service time is optimized. In [START_REF] Khawam | Individual vs. Global Radio Resource Management in a Hybrid Broadband Network[END_REF], authors studied the global approach where the system load balances the downlink traffic of every user between two RATs WiMAX and WiFi in a way to privilege the overall system performances. The hybrid approach [START_REF] Akyildiz | A ubiquitous mobile communication architecture for next-generation heterogeneous wireless systems[END_REF] combines the two previous ones and benefits from feedback in order to enhance local decision. In this paper, we adopt the global approach where the system allocates the downlink traffic between two different technologies. Based on our work in [START_REF] Khawam | Individual vs. Global Radio Resource Management in a Hybrid Broadband Network[END_REF], we formulate this approach as an optimization problem that consists in minimizing the total delay defined as the sum of all individual delays. Then we study the case where we have two types of WiFi [START_REF]Wireless LAN Medium Access Control (MAC) and Physical Layer (PHY) specifications[END_REF] technologies (802.11b and 802.11g) as if we had two different RATs. The challenge we face here is the high computational complexity to obtain the optimal solution for a large number of users. Therefore, as a contribution to our previous work where we just focused on the formulation of the RAT selection problem and gave a solution based on game theory, we propose here four distributed heuristic algorithms for RAT selection. Then we evaluate by simulation the performance of these algorithms with respect to the optimal one. Three different scenarios are considered regarding the spatial users distribution for this evaluation.

The rest of the paper is organized as follows. The system model is described in Section II. The studied 802.11b and 802.11g scenario is described in Section III. The proposed heuristic algorithms are presented in Section IV. Simulation results are given in Section V. Conclusions are given in Section VI.

II. SYSTEM MODEL ON THE DOWNLINK AND ASSUMPTIONS

We consider a geographic region that is entirely covered by a primary RAT and partly covered with several hotspots of a secondary RAT as shown in Figure 1. Only mobile users present in secondary RAT hotspots benefit from load balancing their traffic between the two available RATs. The other mobile users profit only from the primary RAT coverage. Also we assume that we are in a static state of the network and users are stationary. In other words, we take a snapshot of a dynamic system and optimize its current state; hence, we consider that the calculation is performed for each user arrival/departure. Furthermore, we assume that the system is in a saturation state which means that we treat a worst case scenario where every user has persistent traffic.

A. Network Structure

The index x is used throughout the paper to designate a given cell, x = 0 for the primary RAT and x = 1, ..., H for the H secondary RAT hotspots. We denote by n x the number of users able to connect to cell x. Since all users can connect to the primary RAT, n 0 is the total number of mobile users and n = n 0 -H x=1 n x is the number of users that are only serviced by the primary RAT. Thus, only n = H x=1 n x users are engaged in the RAT selection; such users will be indexed by k = 1, ..., n and will be coined hybrid users. We term by h(k) the secondary RAT hotspot of user k. Note that even if we have several hotspots h(k), a user k can be connected only to one h(k).

B. Traffic Model

We denote by θ k the instantaneous fraction of time (or the percentage of traffic) during which user k is assigned to its secondary RAT hotspot. Hence, 1-θ k is the remaining fraction of time spent by user k in the primary RAT. We assume the existence of a central entity responsible of routing the downlink traffic of each user either to the primary or secondary RAT depending on the user association.

C. Cost Function

We denote by T k,x , the cost function of user k when sending all its traffic in RAT x. T k,x will be a function of the percentage of traffic of all other users θ -k in the cell. Thus the cost in the primary RAT is:

T k,0 = f 0 (1 -θ -k ) (1)
and in the secondary RAT is:

T k,h(k) = f h(k) (θ -k ) (2)
Therefore, the cost function of a hybrid user k, denoted by C k (θ k ) is given by:

C k (θ k ) = T k,0 • (1 -θ k ) + T k,h(k) • θ k (3)

D. Non Linear Problem

In the global approach, the system load balances the downlink traffic of each hybrid user among RATs in order to optimize the total cost. Thus, this approach can be formulated as an optimization problem (P) that consists in minimizing the total network cost. This cost, denoted by C t , is defined as the sum of all individual user costs and is given by:

C t (θ k , k ∈ N ) = n k=1 C k (θ k )
, where N denotes the set of hybrid users. Consequently, the optimization problem (P) is given by:

(P) : Minimize C t (θ k , k ∈ N ) = n k=1 C k (θ k ) Subject to: 0 ≤ θ k ≤ 1, k = 1, .., n

III. HETEROGENEOUS WIFI SCENARIO

The model we developed above is applicable to different types of technologies. As the 802.11 WLANs are being widely deployed, we study the case where we have two kinds of WiFi technologies as if we had two different RATs, considering the 802.11g technology as the primary RAT and the 802.11b technology as the secondary one.

A. Radio Conditions

In each cell, the peak rate of each user depends on its radio conditions. Let χ k,x denote the instantaneous peak rate of user k in cell x. Each hybrid user perceives two peak rates χ k,0 and χ k,h(k) depending on the chosen RAT.

B. Data Rate Model in WiFi

As we consider only the downlink traffic and neglect the 802.11 waiting times (i.e DIFS, SIFS), the CSMA/CA access scheme leads to a fair rate sharing of the channel among users. Therefore, all users have the exact same throughput. The throughput in question is less than the minimum peak rate that a single user can perceive when alone in the cell [START_REF] Heusse | Performance Anomaly of 802.11b[END_REF]. Therefore, the data rate of user k in 802.11g is:

R k,0 = 1/ n0 i=1 1 l {user i admitted in 802.11g} χ i,0 (4)
where

1
R k,h(k) = 1/ n i h(i)=h(k) 1 l {user i admitted in 802.11b hotspot h(k)} χ i,h(i) (5)

C. Cost Function

We consider that T k,x , the cost function of user k in cell x adopted in the theoretical fair rate model, is the amount of time necessary to send a data unit in cell x. In other words, the cost function is the inverse image of the throughput. In 802.11g, we deduce from Equation 4 the following:

T k,0 = E[1/R 0 ] = 1 χ k,0 + n0 i=1,i =k 1 χ i,0 + n i=1,i =k 1 -θ i χ i,0 (6)
The same applies in 802.11b hotspot h(k) using Equation 5:

T k,h(k) = E[1/R k,h(k) ] = 1 χ k,h(k) + n i=1,i =k h(i)=h(k) θ i χ i,h(i) (7)
Therefore, the cost function of a hybrid user k given in Equation 3 will be its expected time necessary to send a unit of data in the hybrid environment.

D. Optimal Solution

In the previous section, the global performance is formulated as the non linear optimization problem (P). In [START_REF] Khawam | Individual vs. Global Radio Resource Management in a Hybrid Broadband Network[END_REF], the authors demonstrated that the global optimal solution for this problem is always reached when θ k = 1 or θ k = 0 for each user k, k = 1, ..., n. Therefore each user is assigned either to the 802.11g or to the 802.11b hotspot.

The time complexity for computing the optimal solution is O(2 n ), where is the number of users. Moreover we need a central entity that has an overall knowledge of the network. Therefore, we are going to introduce distributed heuristic algorithms for RAT selection in the next section.

IV. DISTRIBUTED HEURISTIC ALGORITHMS

Due to the high computational complexity of the global optimal solution, we propose four distributed heuristic algorithms for RAT selection starting with a simple algorithm based on the distance between the user and each of the primary and secondary access points (APs). In the second one, in order to give the possibility to distribute users between the two technologies, we add uncertainty to the selection so that the probability to connect to one technology is inversely proportional to the distance between the user and the corresponding AP. The third algorithm is based on the peak rate that each user receives from the two technologies. Finally in the last algorithm, to not connect all users to the broader technology, we combine peak rate and probability in a manner the probability to connect to one technology is proportional to the peak rate received by the user from this technology.

A. Distance Based Selection

In this algorithm we take into consideration the distance between the user and each of the two APs where the user connects to the nearest AP. In fact, this distance reflects the power of the received signal at the user side. Algorithm 1 describes the different steps of the distance based selection algorithm. the user will be connected to 802.11b, 4: else the user will be connected to 802.11g. the user will be connected to 802.11g as the probability to connect to 802.11g is:

d b d b + d g (8)
5: else the user will be connected to 802.11b as the probability to connect to 802.11b is:

d g d b + d g (9)
6: end if Algorithm 3 Peak rate based algorithm 1: Each user measures its received peak rates χ b and χ g respectively obtained from 802.11b and 802.11g APs. 2: if χ b < χ g then 3:

the user will be connected to 802.11g, 4: else the user will be connected to 802.11b. 5: end if

B. Probabilistic Distance Based Selection

As mentioned earlier, in this algorithm we combine the distance to a probability function in a way that the probability to connect to a given technology is inversely proportional to the distance between the user and the corresponding AP. Therefore, the closer the user is to the AP, the higher its probability to be connected to this AP. The different steps considered in this algorithm are shown in Algorithm 2.

C. Peak Rate Based Selection

This algorithm is based on the peak rate each user receives from the two technologies in such a way the user connects to the RAT where it gets the maximum peak rate. The different steps of this algorithm is described in Algorithm 3.

D. Probabilistic Peak Rate Based Selection

In this algorithm, we combine the peak rate with a probability function. We obtain that the probability to connect to Algorithm 4 Probabilistic peak rate based algorithm 1: Each user measures its received peak rates χ b and χ g respectively obtained from 802.11b and 802.11g APs. 2: Each user draws a random variable x following a uniform distribution between 0 and χ b + χ g . 3: if x < χ b then 4: the user will be connected to 802.11b as the probability to connect to 802.11b is:

χ b χ b + χ g ( 10
)
5: else the user will be connected to 802.11g as the probability to connect to 802.11g is:

χ g χ b + χ g ( 11
)
6: end if a given technology is proportional to the peak rate received by the user from the corresponding AP. Thus, the higher the peak rate received from the AP, the higher will be the user probability to be connected to this AP. The different steps considered in this algorithm are shown in Algorithm 4.

In order to evaluate the efficiency of these algorithms we will compare by simulation the global optimal solution to the results obtained by these heuristic algorithms in the next section.

V. SIMULATION RESULTS

Simulations in this work are done using NS2 Network Simulator. We consider a geographic region that is entirely covered by 802.11g and partly covered by a 802.11b hotspot. There are 12 hybrid users generated randomly following a uniform distribution. These users are downloading a Pareto traffic over User Datagram Protocol (UDP) with a Burst Time of 400 ms, an Idle Time of 100 ms, a source Data Rate of 20Mb/s (which leads to saturation state of the network according to the assumption presented in section II) and an average Packet Size of 1000 Bytes. We consider three different scenarios regarding the spatial users distribution with respect to the WiFi APs:

• Scenario I: 50% of users are closer to the 802.11g AP.

• Scenario II: 75% of users are closer to the 802.11g AP.

• Scenario III: 25% of users are closer to the 802.11g AP. For each scenario, five users association schemes are simulated. The first scheme implements the optimal global solution while the four others correspond to the heuristic approaches. More precisely, for each scenario, we calculate the total cost defined as the sum of all individual costs for the optimal solution and compare it with the costs obtained by the proposed heuristic algorithms. For simplicity, we denote by O the simulated global network cost for the optimal solution, by D the distance based solution, by P D the probabilistic distance based solution, by R the peak rate based solution, and by P R the probabilistic peak rate solution. Therefore, we consider 10 ratios that drive our comparison as shown in table (I).

A. Sceanrio I

In this scenario we consider that 50% of users are closer to the 802.11g AP. For all the considered ratios given in table (I), we show the 95% confidence interval for scenario I in Figure 2. The intervals in blue represent the ratios of the proposed Figure 2. Compared costs for scenario I algorithms costs with respect to the optimal one (Ratios 1, 2, 4 and 7). The red intervals represent the ratios of the heuristic costs with respect to the distance based solution (Ratios 3, 5 and 8). The green ones represent the ratios with respect to the probabilistic distance based (Ratio 6 and 9). Finally, the purple interval represents the cost of probabilistic peak rate cost over the peak rate one (Ratio 10). Note that when the ratio interval is higher than 1, we can conclude that the solution corresponding to the denominator is more efficient than the one corresponding to the numerator and vice versa. Simulation results confirm that the optimal solution is the most efficient solution that is why all blue intervals are above 1. Moreover, we see that ratio 7 has an average value of 1.2 so we conclude that the probabilistic peak rate based solution in this scenario is the closest to the optimal one. For the red intervals, ratio 3 interval is above 1 so the distance based solution here is more efficient than the probabilistic distance one since users are in this scenario evenly distributed among the two APs. Ratio 5 has an average value almost close to 1 with a narrow interval. Therefore, we consider that the peak rate based and the distance based solutions present the same efficiency. This result is expected since the 2 APs have the same transmit power, which implies that distance and peak rate are equivalent. From the purple interval (below 1), we affirm that probabilistic peak rate based solution is more efficient than the peak rate based solution. This is due to the fact that the peak rate based solution causes the underutilization of the 802.11b technology while the probabilistic peak rate based solution avoids the rush to the AP that offers better peak rate by adding uncertainty to the selection so that not all users will be connected to the broader technology.

B. Sceanrio II

In scenario II, we proceed as in scenario I but here we consider different spatial user distribution so that 75% of users are closer to the 802.11g AP. The obtained results are shown in Figure 3. The blue intervals are all above 1 and confirm that Figure 3. Compared costs for scenario II the optimal solution is again the most efficient one. Moreover, we see that ratio 1 has an average value close to 1 therefore the distance based solution gives results closest to the optimal one in this scenario. This is predictable since 75% of users are closer to the 802.11g AP, thus they will be well served by this AP given the high capacity offered by 802.11g technology. For the green intervals, they are both under the value 1, therefore in this scenario the peak rate based and the probabilistic peak rate based solutions outperform the probabilistic distance one.

C. Sceanrio III

In this scenario, we consider that 75% of users are closer to the 802.11b AP. We obtain the results shown in Figure 4. Ratio 1 has an average value of 2.36 therefore the distance based solution in this scenario is not efficient compared to the optimal one. In fact, with 75% of users connected to 802.11b AP, congestion will occur causing degradation of users QoS. Ratio 2 has an average value of 1.4 with a narrow confidence interval, so we can say that the probabilistic distance based solution gives us here good results. Consequently, by combining probability with distance in this scenario we give the possibility to users to connect to the 802.11g technology and therefore enhance global performance. For ratio 9, the entire confidence interval is under the value 1, thus the probabilistic peak rate based solution for this scenario is more efficient than the probabilistic distance one.

VI. CONCLUSION

In this paper, we described a hybrid system model for RAT selection where the system manages the users RAT association to enhance global performance. Then we studied the case where we have two types of WiFi technology (802.11b and 802.11g) as if we had two different RATs, to find that in the optimal solution of this approach each user is connected to a single RAT. Due to the high computational complexity of the global optimal solution for a large number of users, we proposed four distributed heuristic algorithms for RAT selection. We concluded that each heuristic algorithm presents efficient results depending on the users distribution so it is interesting to implement an entity that chooses the appropriate algorithm. Mainly, the probabilistic peak rate based solution gave us results close to the optimal one. Moreover, a comparison between the proposed algorithms was provided in the three different scenarios. For future work, we suggest to implement these algorithms to the WiMAX -WiFi scenario since WiMAX offers high data rates. This will be an extension of the current work because the model we developed for the global approach treat the fair rate scheduler model for the primary and secondary RAT, thus it is interesting to evaluate these heuristics in the case of the fair time or fair rate -fair time hybrid schedulers since the global optimal solution will boil down to a unique RAT selection.

Figure 1 .

 1 Figure 1. Hybrid wireless access environment

Algorithm 1

 1 Distance based algorithm 1: Each user measures the distances d b and d g separating it respectively from the 802.11b and 802.11g APs. 2: if d b < d g then 3:

Figure 4 .

 4 Figure 4. Compared costs for scenario III

  Each user measures the distances d b and d g that separate it respectively from the 802.11b and 802.11g APs. 2: Each user draws a random variable x following a uniform distribution between 0 and d b + d g . 3: if x < d b then

	5: end if
	Algorithm 2 Probabilistic distance based algorithm
	1: 4:


Ying Hu

Mi Youn Jung

Ahmed Oukili

Guanyu Yang

Jean-Claude Nunes
email: jean-claude.nunes@univ-rennes1.fr

Jérôme Fehrenbach

Gabriel Peyré

Marc Bedossa

Limin Luo

J.-C Nunes

C Toumoulin

Fellow, IEEE. L D Cohen

J C Nunes

Sparse reconstruction from a limited projection number of the coronary artery tree in X-ray rotational imaging

Keywords: reconstruction, Maximum a posteriori (MAP), minimal path, non-local active contours, X-ray rotational coronary angiography

à 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.

I. INTRODUCTION

Rotational X-ray imaging devices provide a new way to explore the possibly existing artery pathologies enhanced after injection of a dye product, the complete range of projections (80 up to 150), giving a pseudo 3-D view of the coronary tree with all the structures moving all together. However the challenge remains today the 3D coronary reconstruction that will provide the ground for a platform dedicated to the planning and execution of percutaneous coronary interventions. Due to the rotation of the source detector system, all the structures are moving and their features (including the most basic one, the intensity, e.g the attenuation) are varying over time and space. The object of interest, the coronary tree, is patient dependent (changes in its branching structure, complex-shaped vessels), has low nonstationary contrast and crossings, superimpositions in the image sequence planes represent additional difficulties to deal with. The coronary motions, determined by the myocardium contraction, can not be reduced to rigid transformations: strong deformations occur with slow and fast phases and movement inversions.

Different approaches have been proposed for this reconstruction depending on whether we consider the set of available projections or only the one acquired at the same phase of the cardiac cycle [START_REF] Schoonenberg | Three-Dimensional coronary visualization, Part 2: 3D reconstruction[END_REF]. In the first option, the 3D+T motion is first estimated throughout the cardiac cycle to perform a motion compensated tomographic reconstruction at a given instant of the cardiac cycle using all the available projections [START_REF] Schäfer | Motion compensated and gated cone beam filtered back projection for 3D rotational X-ray angiography[END_REF]. The construction of this motion model generally relies on modeling techniques to extend the information to the time domain and is most often estimated from extracted 2D centerlines [START_REF] Bousse | Motion compensated tomography reconstruction of coronary arteries in rotational angiography[END_REF]. In the second option, the reconstruction is carried out from a limited number of views (typically 4 -6 projections), corresponding to a single cardiac phase selected by ECG-gating. This reconstruction can be obtained by solving a static tomographic inverse problem. But the low number of projections leading to a severe angular undersampling, prior information is often introduced to improve the quality and accuracy of this reconstruction [START_REF] Hansis | Evaluation of iterative sparse object reconstruction from few projections for rotational coronary angiography[END_REF]. Other methods still make use of the epipolar geometry constraints and feature matching techniques to find corresponding structures in each view and recover their geometry [START_REF] Zheng | Sequential reconstruction of vessel skeletons from X-ray coronary angiographic sequences[END_REF].

We address here this illness problem considering background subtracted projections so that a vessel only 3D reconstruction can be performed, introducing thus a sparse representation of the data in the space domain. The reconstruction procedure goes through a preliminary segmentation stage and a l p norm regularization (p < 2) within a statistical MAP reconstruction algorithm. These two methods are described in Section II. Preliminary results are given in Section III on simulated data.

II. METHOD

The segmentation algorithm takes benefit of both minimal path and level set evolution methods to extract the coronary tree on each projection of the sequence. Then, the background is estimated to build background subtracted images and deal with sparse data. We address respectively, for comparison, the l 0 and l 1 norm regularizations within a maximum a posterior (MAP) estimation. The optimization process relies on a separable paraboloidal surrogate (SPS) algorithm which was introduced by Erdogan et al. for PET image reconstruction under Poisson statistics [START_REF] Erdogan | Ordered subsets algorithms for transmission tomography[END_REF] and later for MSCT reconstruction [START_REF] Hu | L0 constrained sparse reconstruction for multi-slice helical CT reconstruction[END_REF]. SPS functions are applied here in the frame of a gaussian distribution data model.

The minimal path technique introduced by Cohen and Kimmel [START_REF] Cohen | Global minimum for active contour models: A minimal path approach[END_REF] has been extensively used for the extraction of tree or tubular tree structures. This approach has several advantages such as finding global minima, fast computation, easy implementation, and more powerful incorporation of user input. On the other hand, level set evolution techniques have a nice advantage, which is that they capture and represent the boundary of object directly, not just some path running through their interior as in the minimal path approach. Recently, an active contour model [START_REF] Jung | Nonlocal active contours[END_REF] was proposed in the level set framework, which constrains the local homogeneity of image features. This is crucial to overcome the difficulty of segmentation due to image inhomogeneity that is often seen in medical images.

1) Minimal Paths: The minimal path technique characterizes a boundary extraction approach that globally minimizes a geodesic active contour energy between two user supplied end points. The model is formulated as

min γ∈Ax 0 ,x E(γ) = L 0 w + P (γ(s))ds = L 0 P (γ(s))ds
where A x0,x is the set of all paths linking x 0 and x, s represents the arc-length, and w > 0 is a constant imposing regularity on the curve. P > 0 is a potential cost function computed from the image, which takes lower values on the interesting features of the image.

To compute the solution associated with the source point x 0 of this problem, a Hamiltonian approach was proposed in [START_REF] Cohen | Global minimum for active contour models: A minimal path approach[END_REF]. This finds the geodesic weighted distance U 0 by solving the Eikonal equation: ∇U 0 (x) = P (x), ∀x ∈ Ω, with U 0 (x 0 ) = 0. Then, the minimal path γ can be retrieved with a simple gradient descent on U 0 from x to x 0 , by solving the ordinary differential equation: dγ(s) ds = -∇U 0 (γ(s)) with γ(0) = x. To solve the Eikonal equation, one can use the Fast Marching algorithm [START_REF] Cohen | Global minimum for active contour models: A minimal path approach[END_REF] due to its lower complexity.

2) Non-local Active Contours: The non-local active contours model [START_REF] Jung | Nonlocal active contours[END_REF] aims at finding a contour that represents the boundary of object of interest, and minimizes the energy incorporated with a level set function φ : Ω → R, , a patch p x around the pixel x and a metric d(•, •) ≥ 0 that accounts for the similarity between patches, and where

F (φ) = Ω Ω ρ(H(φ(x)), H(φ(y)))K(x, y)dxdy + γL(φ) where K(x, y) = G σ (x -y)d(p x , p y ) with G a (t) = e -t 2 /2a
L(φ) = Ω ∇H(φ(x)) dx. H is a smoothed Heaviside function, and ρ is an indicator function such that ρ(u, v) = 1 if u = v, 0 otherwise (e.g. ρ(u, v) = uv+(1-u)(1-v)).
The parameter σ > 0 controls the scale of the local homogeneity one requires for the segmented object. The second term enforces the smoothness of the contour, thus γ > 0 is a regularization parameter.

The solution of the minimization problem is computed by solving the evolution equation for an artificial time t ≥ 0:

∂φ(t,x) ∂t = -∇F (φ(t, x)) with φ(0, x) = φ 0 (x)
. See [START_REF] Jung | Nonlocal active contours[END_REF] for more details.

3) Proposed Segmentation Algorithm: In the proposed algorithm, the minimal path technique is first used for an initial estimation of the tubular tree structure. The non-local active contours model is then performed, to make a level set function constructed from the minimal paths move towards the boundary. Thus, the segmentation algorithm is as follows:

• (Estimation step) Given x 0 and end points {x i } k i=1 , • Compute minimal paths γ x0,xi (i = 1, ..., k).

• Compute an initial level set function φ 0 as a signed distance function along the minimal paths. • (Correction step) Solve the non-local segmentation model starting with the initial level set function φ(0, •) = φ 0 .

B. 3D reconstruction stage

Rotational imaging provides a sequence of 2D projections acquired under continuous rotation of the C-arm. The acquisition trajectory is limited to a circular arc with a fixed caudo-cranial angulation. Let consider Y = {Y 1 , Y 2 , ..., Y L } to be the projection images acquired at primary angles θ = {θ 1 , θ 2 , ..., θ L }, L being the number of considered rotation angles. If we assume the cardiac motion is relatively regular and periodical, then through ECG gating, we can choose a set of projections Y s , which corresponds to the same 3D heart motion (i.e the same cardiac phase s). The 3D coronary tree f s at a given instant s of the cardiac cycle, will be reconstructed from a very few number of projections, which renders the problem strongly ill-posed. We will further simplify the notation by using Y instead of Y s for conciseness. Projections Y are non-subtracted i.e they include both the background tissues and the coronary tree enhanced with the contrast agent. Thus, the expectation of the measurement E(Y ) can be written under the following form:

Ȳi = Ȳbi exp (-[A.f ] i ) (1)
with [Af ] i = j a ij f j (j = 1, ..., J, i = 1, ..., I) and (J, I) being the voxel and pixel numbers in the volume and on the detector plane respectively. A is the cone beam projection operator, a ij denotes the contribution of voxel j in the computation of pixel i, and Y bi is the i th element of the background tissues Y b . The logarithm application on eq. 1 provides a new image that only contains the structure of interest: ḡi = [A.f ] i with ḡi = -log Ȳi Ȳb i . The data are thus supposed to follow a Gaussian distribution model:

P (g i |f ) = 1 σ i √ 2π exp - 1 2
g i -g i σ i 2 (2)
1) MAP model: According to the Bayesian theory, the estimate of the unknown object f is computed as a function which maximizes the posterior density P (f |g) and is given by:

f (g) = arg min f ≥0 (-log P (g|f ) -log P (f |β)) (3)
The log-likelihood function log P (g|f ) = L(g|f ) is:

L(g|f ) = - 1 2 (g -ĝ) T Σ -1 (g -ĝ) T -const ( 4
)
where ĝ = Af and Σ is the covariance matrix: Σ = diag{σ i } 2 with σ 2 i = max(σ 2 min , log (ĝ)). The prior function has the form of a Gibbs distribution

P (f |β) = exp (-βR(f )) Z(β)
, Z(β) being a scaling constant. Thus, log P (f |β) = -βR(f ) with:

R(f ) = J j=1 k∈Nj ω jk ψ(f j -f k ) ( 5
)
where N j defines the neighborhood of the j th voxel, and ω jk is a positive value that expresses the interaction degree in clique kj. ψ(t) is a symmetric function that penalizes the pairwise differences between neighbouring voxels. Various kinds of potential function ψ have been proposed in the literature that are l 2 norm [START_REF] Hansis | Evaluation of iterative sparse object reconstruction from few projections for rotational coronary angiography[END_REF], l 1 norm [START_REF] Bousse | Motion compensated tomography reconstruction of coronary arteries in rotational angiography[END_REF] and l 0 norm [START_REF] Hansis | Evaluation of iterative sparse object reconstruction from few projections for rotational coronary angiography[END_REF]. We compare 2 of these priors: l 0 norm and l 1 norm respectively defined by:

ψ(t) = |sgn(t)| and ψ(t) = |t|.
The l 0 norm is non convex and not continuous around zero. This brings some difficulty for the optimization. A way for solving it is to go through surrogate Functions. Any function satisfying: lim ρ→0 ψ(t, ρ) = |sgn(t)| can be taken as surrogates, such as ψ(t, ρ) = 2 π arctan |t| ρ .

2) Optimization:

The key of the SPS method is to find a decomposition of the objective function into a simple quadratic form that simplifies its optimization, guaranties a faster global convergence and reduces the computational time [START_REF] Zheng | Sequential reconstruction of vessel skeletons from X-ray coronary angiographic sequences[END_REF]. We define:

L(ĝ) = i L i ( ĝi ) = i (gi-ĝi) 2 2σi 2 . Then Lj (f j , f ) = i L i (g i , ĝi
). The voxel based quadratic likelihood function is then given by Qj (

f j ) = Lj ( f ) + Lj (f j )(f j -fj )-1 2 d j (f j -fj ) 2 with Lj (f j ) = i a ij Li (ĝ i ) and d j = i a ij 2 σ i 2 .
In the case of the non convex prior, we apply an update into 2 steps:

f j k+1,1 = max(0, f j k - Qk j Qk j )
and f j k+1,2 = max(0, f j k+1,1 + β ∂R(f ) ∂fj ).

III. RESULTS

Evaluations have been carried out on simulated data obtained from a 3D dynamic sequence acquired on a 64-slice GE LightSpeed CT scan [START_REF] Yang | A multiscale tracking algorithm for the coronary extraction in MSCT angiography[END_REF]. A sequence of 20 3D binary coronary trees was built corresponding to 20 different cardiac phases. A C-arm rotational R-X coronary angiography was then simulated using the Siemens Axiom System imaging protocol. The detector plane (200mm) 2 was uniformly sampled into 512 pixels. Reconstructions were performed within a volume of (110mm) 3 , 80 projections of the 3D binary tree were generated, uniformly spaced over the range RAO 90 • to LAO 90 • . 4 cardiac cycles were considered, this means that a volume (associated with a phase s of the cardiac cycle) was projected 4 times during the acquisition, according to different viewpoints. The projection operator A has been computed according to [START_REF] Man | Distance-driven projection and backprojection in three dimensions[END_REF]. We used then the method of low order polynomials approximation to build background tissues images from the CT scan datasets. Fig. 1 displays four simulated projections Y s , s being the phase 0 and the projection angles RA0 90 • , RAO 30 • , LAO 30 • and LAO 90 • .

A. Segmentation results

The fig. 2

B. Reconstruction results

Background images Y b were then estimated by subtracting the segmented images (fig. 3) from original projections Y and applying an inpainting method on Y b to interpolate the missing background data (at the coronary tree location). The background segmented images g are then obtained by the logarithmic subtraction (log Y b -log Y ).

The hyperparameter β = 0.002 for l 0 norm, β = 0.0005 and for l 1 norm. We applied a relaxation scheme to estimate the optimal value ρ in the l 0 norm prior computation. Fig. 4 displays the 3D coronary reconstruction obtained from background segmented images. The reconstruction quality is evaluated using DICE overlap criterion. This criterion is equal to 0.42 for both l 0 and l 1 norms. For comparison, we performed the reconstruction of ideal coronary trees directly obtained from the 3D reference sequence and got a DICE of 0.99 for l 0 norm and 0.97 for l 1 norm. This result shows that the reconstruction algorithm is highly dependent on the segmentation-based background estimation.

IV. CONCLUSION

The coronary artery tree reconstruction from a very few number of projections appears to be a very strong ill-posed problem. Indeed the cardiac motion leads to a severe angular undersampling. Moreover, we restricted the projection number to four for the static reconstruction in order to position ourself in situation of clinical routine where the acquisition is made over a period of 4 seconds. We proposed to solve this problem by simplifying the image contents with a preliminary segmentation stage, which is based on minimal paths and non local active contours to build background subtracted projections. We compared then l 0 and l 1 norm regularizations within a statistical MAP reconstruction algorithm to deal with the spatial sparsity of the data. Evaluations have been performed by means of a realistic phantom representing an arterial tree extracted from a sequence of MDCT datasets in order to get a ground truth. Preliminary results make appear that the reconstruction algorithm is very sensitive to the quality of the segmentation and the background estimation. Our further objective will thus be to improve the background estimation as well as to work on a new prior that will better discriminate the coronary/background structures during the reconstruction process.

Fig. 1 .

 1 Fig. 1. Simulated projections at phase 0. Projection angles are respectively (RA0 90 • , RAO 30 • , LAO 30 • and LAO 90 • ).

Fig. 2 .

 2 Fig. 2. Segmentation results. (1st) original image, (2nd-4th) minimal paths, resulting curves of NLAC, extracted vessels. Blue and green points represent a starting point and end points respectively.

  presents segmentation results of the algorithm described in II-A. For the potential P , we computed an enhanced image f 0 = G b * f -f with b = 10, which cleaned the original image f by subtracting the smoothly varying background. We used P (x) = w + |f 0 (x)max x f 0 (x)| with w = 0.001. For the non-local active contours, we used the L 2 distance d based on the image f with patches of size 3 × 3 pixels, G σ of size 15 × 15 pixels with σ = 1000, and we set γ = 40. The curve evolution method could correct erroneously segmented (or missing) parts that are initially estimated by the minimal path technique.

Fig. 3 .

 3 Fig. 3. Segmented images g for the four projections of the phase 0 of the cardiac cycle.

Fig. 4 .

 4 Fig. 4. Coronary tree reconstruction from the 4 projections depicted in the fig. 1. First line: reconstruction from the ideal tree (left: with l 0 norm, right: with l 1 norm). Second line: reconstruction from the background segmented images (left: with l 0 norm, right: with l 1 norm).

  2

V. ACKNOWLEDGMENTS This work was supported by ANR grant MESANGE ANR-08-BLAN-0198 and the National Natural Science Foundation NO. 31100713.