Thierry Denoeux

Nicole El Zoghby

Véronique Cherfaoui

Antoine Jouglet

Optimal object association in the Dempster-Shafer framework

Keywords: Belief functions, Evidence theory, data fusion, assignment problem

come

I. INTRODUCTION

Object association refers to the task of matching two finite sets of objects E = {e 1 , . . . , e n } and F = {f 1 , . . . , f p }, with possibly different cardinalities. This problem arises, for instance, in multiple target tracking applications [START_REF] Bar-Shalom | Tracking and data association[END_REF], [START_REF] Blackman | Design and Analysis of Modern Tracking Systems[END_REF], in which we need to estimate the status of mobile objects (such as targets or storm cells [START_REF] Einfalt | A radar rainfall forecasting method designed for hydrological purposes[END_REF]) that are detected at different times by a single sensor. Data association then consists in deciding which observation should be used to update each track. In this case, the two sets of objects are a set of tracks and a set of observations. Another class of problems in which object association is needed is sensor fusion (see, e.g., [START_REF] Houenou | A Track-To-Track Association Method for Automotive Perception Systems[END_REF]). In that case, objects are typically perceived by two sensors and we need to perform observation-to-observation or track-to-track association, depending on the level of sensor outputs. The study of the multi-dimensional case, in which we have more than two sets of objects (perceived, e.g., by several sensors), is deferred until the end of this paper.

Usually, object association is performed under the assumption that each object in one set should be matched with at most one object in the other set. An object in E may have no counterpart in F because it has disappeared between two successive time frames, or because it has not been perceived by one of the sensors. Mathematically, we are thus searching for a relation R ⊆ E × F such that, for all i, j and k:

(e i , f j ) ∈ R and (e i , f k ) ∈ R ⇒ j = k (1a)
and

(e i , f k ) ∈ R and (e j , f k ) ∈ R ⇒ i = j. (1b)
Any such relation may be described by a matrix R of size (n, p) such that R ij = 1 if (e i , f j ) ∈ R and R ij = 0 otherwise (by abuse of notation but without any risk of confusion, we use the same notation for the relation and its corresponding matrix).

In this paper, this problem will be investigated in the framework of the theory of belief functions, also referred to as Dempster-Shafer theory or Evidence theory [START_REF] Shafer | A mathematical theory of evidence[END_REF]. We assume that we receive evidence about the possible association of each object pair (e i , f j ). Mathematically, such evidence can be represented by a mass function m ij on the frame Θ ij = {0, 1}, such that m ij ({1}) = α ij is the probability of knowing that R ij = 1, m ij ({0}) = β ij is the probability of knowing that R ij = 0, and m ij ({0, 1}) = 1 -α ij -β ij is the probability of knowing nothing at all about R ij .

Based on such evidence, we would like to choose the "best" relation R * , among the set R of all relations verifying (1a)-(1b). The problem of selecting a relation in R based on pairwise mass functions m ij has been addressed by several authors (see, e.g., [START_REF] Rombaut | Decision making in data fusion using Dempster-Shafer's theory[END_REF]- [START_REF] Dallil | Sensor fusion and target tracking using evidential data association[END_REF]). However, only heuristic solutions have been provided until now. The most elaborate solution so far, proposed by Mercier et al. [START_REF] Mercier | Object association with belief functions, an application with vehicles[END_REF], consists in first combining the mass functions {m ij } p j=1 for each i, and then finding the relation R with maximum pignistic probability [START_REF] Smets | The Transferable Belief Model[END_REF]. However, this algorithm involves enumerating all the elements of R, which quickly becomes intractable when n and p are not very small. Additionally, the method lacks a fundamental symmetry property, as it may give different results if the sets E and F are interchanged.

In this paper, we show that the above problem, formalized as the search for the most plausible relation R * in R, can be transformed into an equivalent linear assignment problem and solved exactly in polynomial time. The method is studied experimentally using both simulated and real multi-sensor data from an intelligent vehicle application. Finally, we show that the three-dimensional extension of this problem, with three object sets E, F and G, is N P-hard.

The rest of this paper is organized as follows. After recalling the necessary background on the theory of belief functions in Section II, the object assignment problem is formalized and solved in Section III. The sensor fusion application is then addressed in Section IV, where experimental results with simulated and real data are also reported. Finally, the three-dimensional version of the object association problem is studied in Section V and Section VI concludes the paper.

II. BACKGROUND ON BELIEF FUNCTIONS

The theory of belief functions has two main components: equivalent representations of a body of evidence (in the form of mass, belief and plausibility functions), and a combination rule for combining independent items of evidence. These two components are reviewed in Subsection II-A and II-B, respectively.

A. Representation of evidence

The theory of belief functions is a framework for reasoning under uncertainty based on the modeling of evidence [START_REF] Shafer | A mathematical theory of evidence[END_REF]. More precisely, let us assume that we are interested in the value of some variable θ taking values in a finite domain Θ, called the frame of discernment. Uncertain evidence about θ may be represented by a (normalized) mass function m on Θ, defined as a function from the powerset of Θ, denoted as 2 Θ , to the interval [0, 1], such that m(∅) = 0 and A⊆Θ m(A) = 1.

(

) 2
Each number m(A) is interpreted as a degree of belief attached to the proposition θ ∈ A and to no more specific proposition.

As argued by Shafer [START_REF] Shafer | Constructive probability[END_REF], the meaning of such degrees of belief can be better understood by assuming that we have compared our evidence to a canonical chance set-up. The setup proposed by Shafer consists of an encoded message and a set of codes Ω = {ω 1 , . . . , ω n }, exactly one of which is selected at random. We know the list of codes as well as the chance p i of each code ω i being selected. Decoding the encoded message using code ω i produces a message of the form "θ ∈ A i " for some A i ⊆ Θ. Then

m(A) = {1≤i≤n:Ai=A} p i (3)
is the chance that the original message was "θ ∈ A". Stated differently, it is the probability of knowing that θ ∈ A. In particular, m(Θ) is, in this setting, the probability that the original message was vacuous, i.e., the probability of knowing nothing.

The above setting thus consists of a set Ω, a probability measure P on Ω and a multi-valued mapping Γ : Ω → 2 Θ \{∅} such that A i = Γ(ω i ) for each ω i ∈ Ω. This is the framework initially considered by Dempster in [START_REF] Dempster | Upper and lower probabilities induced by a multivalued mapping[END_REF]. The triple (Ω, P, Γ) formally defines a finite random set [START_REF] Nguyen | An Introduction to Random Sets[END_REF]. Each piece of evidence can thus be represented by a random set, which induces a mass function.

To each normalized mass function m, we may associate belief and plausibility functions from 2 Θ to [0, 1] defined as follows:

Bel(A) = P ({ω ∈ Ω|Γ(ω) ⊆ A}) = B⊆A m(B) (4a) P l(A) = P ({ω ∈ Ω|Γ(ω) ∩ A = ∅}) = B∩A =∅ m(B), (4b)
for all A ⊆ Θ. These two functions are linked by the relation P l(A) = 1 -Bel(A), for all A ⊆ Θ. Each quantity Bel(A) may be interpreted as the degree to which the evidence supports A, while P l(A) can be interpreted as the degree to which the evidence is not contradictory to A. The following inequalities always hold: Bel(A) ≤ P l(A), for all A ⊆ Θ.

The function pl : Θ → [0, 1] such that pl(θ) = P l({θ}) for all θ ∈ Θ is called the contour function associated to m.

B. Combination of evidence

A key idea in Dempster-Shafer theory is that beliefs are elaborated by aggregating different items of evidence. The basic mechanism for evidence combination is Dempster's rule of combination, which can be naturally derived using the random code metaphor as follows. Let m 1 and m 2 be two mass functions induced by triples (Ω 1 , P 1 , Γ 1 ) and (Ω 2 , P 2 , Γ 2 ) interpreted under the random code framework as before. Let us further assume that the codes are selected independently. For any two codes ω 1 ∈ Ω 1 and ω 2 ∈ Ω 2 , the probability that they both are selected is then

P 1 ({ω 1 })P 2 ({ω 2 }), in which case we can conclude that θ ∈ Γ 1 (ω 1 ) ∩ Γ 2 (ω 2 ). If Γ 1 (ω 1 ) ∩ Γ 2 (ω 2 ) = ∅,
we know that the pair of codes (ω 1 , ω 2 ) could not have been selected: consequently, the joint probability distribution on Ω 1 × Ω 2 must be conditioned, eliminating such pairs [START_REF] Shafer | Constructive probability[END_REF]. This line of reasoning leads to the following combination rule, referred to as Dempster's rule [START_REF] Shafer | A mathematical theory of evidence[END_REF]:

(m 1 ⊕ m 2 )(A) = 1 1 -κ B∩C=A m 1 (B)m 2 (C) (5)
for all A ⊆ Θ, A = ∅ and (m 1 ⊕ m 2 )(∅) = 0, where

κ = B∩C=∅ m 1 (B)m 2 (C) (6)
is the degree of conflict between m 1 and m 2 . If κ = 1, there is a logical contradiction between the two pieces of evidence and they cannot be combined. Dempster's rule is commutative, associative, and it admits as neutral element the vacuous mass function defined as m(Θ) = 1.

Dempster's rule can be easily expressed in terms of contour functions: if pl 1 and pl 2 are the contour functions of two mass functions m 1 and m 2 , then the contour function of m 1 ⊕ m 2 is, using the same symbol ⊕ as used for mass functions and contour functions:

(pl 1 ⊕ pl 2 )(θ) = pl 1 (θ)pl 2 (θ) 1 -κ , (7)
for all θ ∈ Θ.

III. OBJECT ASSOCIATION

In this section, we first show that the problem of finding the most plausible matching between two sets of objects, based on independent evidence pertaining to each pair of objects, can be formalized as a binary linear programming problem (Subsection III-A). We then show in Subsection III-B that this problem is equivalent to a linear assignment problem and we give the complexity of this problem.

A. Problem formalization

As explained in Section I, we assume that the available evidence about the association between the sets E and F consists in np mass functions m ij , 1 ≤ i ≤ n, 1 ≤ j ≤ p. Each m ij encodes a piece of evidence about a binary variable R ij that equals 1 if e i and f j correspond to the same entity, and 0 otherwise. We note that the absence of information about the association between e i and f j may be encoded by the vacuous mass function such that m ij ({0, 1}) = 1. Typically, m ij is based on a measure of similarity between some attributes describing the objects. Concrete examples will be discussed in Section IV.

The key idea behind our approach is to express all the available evidence in the frame of discernment R defined as the set of all possible matchings between E and F verifying (1a)-(1b). Assuming independence between the np items of evidence, the np mass functions can then be combined using Dempster's rule [START_REF] Shafer | A mathematical theory of evidence[END_REF], and the plausibility of any relation R ∈ R may be simply calculated using [START_REF] Gruyer | Matching and decision for vehicle tracking in road situation[END_REF].

Let R ij denote the set of relations that match objects e i and f j :

R ij = {R ∈ R|R ij = 1}. (8)
Each mass function

m ij on Θ ij = {0, 1} may be expressed in R by transferring the mass m ij ({1}) = α ij to R ij , m ij ({0}) = β ij to R ij and m ij ({0, 1}) = 1 -α ij -β ij to R, where R ij denotes the complement of R ij .
Let pl ij denote the corresponding contour function. It has the following expression:

pl ij (R) = 1 -β ij if R ∈ R ij , 1 -α ij otherwise, (9)
for all R ∈ R, which can be expressed more concisely as follows:

pl ij (R) = (1 -β ij ) Rij (1 -α ij ) 1-Rij . ( 10
)
Let m denote the mass function on R obtained by combining the np mass using Dempster's rule. From [START_REF] Gruyer | Matching and decision for vehicle tracking in road situation[END_REF], its contour function pl is proportional to the product of the np mass functions pl ij :

pl(R) ∝ i,j (1 -β ij ) Rij (1 -α ij ) 1-Rij , (11)
and its logarithm is

ln pl(R) = i,j [R ij ln(1 -β ij )+ (1 -R ij ) ln(1 -α ij )] + C, ( 12
)
where C is a constant and it is assumed that β ij < 1 and α ij < 1 for all i and j. We note that the situation where β ij = 1 or α ij = 1 (which corresponds to the case where we have absolute certainty that objects e i and f j should, or should not be matched, respectively) can be easily accounted for by setting β ij = 1or α ij = 1for some arbitrary small > 0.

The most plausible relation R * can thus be found by solving the following linear optimization problem:

max R   i,j w ij R ij   under the constraints (1a)-(1b), with w ij = ln 1 -β ij 1 -α ij . ( 13
)
Before studying this problem in the next section, we may observe that prior knowledge about the true object association can be easily incorporated in this framework. For instance, assume that relations that match more objects are considered a priori as more plausible or, on the contrary, less plausible. Such prior knowledge can be represented by a belief function with contour function pl 0 such that:

pl 0 (R) ∝ exp   λ i,j R ij   , ( 14
)
where λ is a scalar parameter. A positive (respectively, negative) value of λ favors relations R with higher (respectively, lower) cardinality. Combining this contour function with the pairwise contour functions pl ij using Dempster's rule and taking the logarithm, we get

ln pl(R) ∝ i,j [R ij ln(1 -β ij )+ (1 -R ij ) ln(1 -α ij ) + λR ij ] + C. (15)
The obtained maximization problem would then have the same form as above, with w ij now defined as

w ij = λ + ln 1 -β ij 1 -α ij . (16)

B. Problem resolution and complexity analysis

The problem previously formalized and denoted as P hereafter, can be stated as the following integer linear program:

max   n i=1 p j=1 w ij R ij   (17)
subject to

p j=1 R ij ≤ 1 ∀i ∈ {1, . . . , n} (18a)
n i=1 R ij ≤ 1 ∀j ∈ {1, . . . , p} (18b)
R ij ∈ {0, 1} ∀i ∈ {1, . . . , n}, ∀j ∈ {1, . . . , p}, (18c)
where constraints (18a) and (18b) are related to Equations (1a) and (1b), respectively. Without loss of generality, it is considered that n ≥ p (if it is not the case, E and F are interchanged). To be solved, P is reduced to the following problem P :

max n i=1 n j=1 w ij R ij (19) subject to n j=1 R ij = 1 ∀i ∈ {1, . . . , n} (20a)
n i=1 R ij = 1 ∀j ∈ {1, . . . , n} (20b)
R ij ∈ {0, 1} ∀i, j ∈ {1, . . . , n}, (20c)
where ∀i ∈ {1, . . . , n}, ∀j ∈ {1, . . . , p}, w ij = max(0, w ij ) and ∀i ∈ {1, . . . , n}, ∀j ∈ {p + 1, . . . , n}, w ij = 0.

Proposition 1. Let R be an optimal solution of problem P and let R be an optimal solution of problem P , then

n i=1 n j=1 w ij R ij = n i=1 p j=1 w ij R ij . Proof: Suppose that n i=1 n j=1 w ij R ij < n i=1 p j=1 w ij R ij .
A new solution R for problem P is extracted from solution R in the following way. First, R is initialized with

R ij = 0, ∀(i, j) ∈ {1, . . . , n} 2 .
For each i ∈ {1, . . . , n} such that

p j=1 R ij = 1, R ij is set to R ij for all j ∈ {1, . . . , p}. At this point of the building of solution R , note that n i=1 n j=1 w ij R ij = n i=1 p j=1 w ij R ij .
Note also that constraints (20a) and (20b) do not necessarily hold since an object of E has not to be necessarily associated with an object of F in R . In this case, there is exactly the same number of indices i ∈ {1, . . . , n} such that n j=1 R ij = 0 as the number of indices j ∈ {1, . . . , n} such that n i=1 R ij = 0. Then, each i ∈ {1, . . . , n} such that n j=1 R ij = 0 is considered iteratively. We search for the smallest indice j with n k=1 R kj = 0 (knowing that this indice necessarily exists) and R ij is set to 1. If j > p, then w ij = 0 by definition of problem P . If j ≤ p, then necessarily w ij = 0. Indeed, w ij > 0 ⇒ w ij > 0 and by setting R ij to 1 (noticing that this new matching can be added to R since p k=1 R ik = n k=1 R kj = 0), we obtain a new solution for problem P with a higher cost, contradicting that R was optimal. Thus, at the end of each iteration

n i=1 p j=1 w ij R ij = n i=1 p j=1 w ij R ij still holds. Finally, R is also a solution of P with n i=1 p j=1 w ij R ij = n i=1 p j=1 w ij R ij > n i=1 p j=1 w ij R ij contradicting the fact that R is an optimal solution of P . Now, suppose that n i=1 n j=1 w ij R ij > n i=1 p j=1 w ij R ij .
A new solution R for problem P is extracted from solution R in the following way. For all i ∈ {1, . . . , n} and for all j ∈ {1, . . . , p} R ij is set to R ij if w ij ≥ 0 and to 0 otherwise. Note that R is a solution of P since if constraints (20a) and (20b) hold, then constraints (18a) and (18b) also hold. Since, by definition,

w ij < 0 ⇒ w ij = 0 then n i=1 p j=1 w ij R ij = n i=1 n j=1 w ij R ij > n i=1 p j=1 w ij R ij ,
contradicting the fact that R is an optimal solution of P . Now, it is clear that an optimal solution for problem P can be obtained by solving problem P . Indeed, from any optimal solution R of problem P , we can build an optimal solution R for problem P with the same cost value (and then optimal) with R ij = R ij if w ij > 0 and R ij = 0 otherwise. Note that building problem P runs in O(n 2 ) time and that building an optimal solution R of problem P from an optimal solution of problem P runs in O(np). It now remains to know how to efficiently solve problem P .

Fortunately, in the operations research literature, problem P is known as the assignment problem and can be solved in O(n 3 ) times with the Hungarian Method [START_REF] Kuhn | The Hungarian method for the assignment problem[END_REF] or in O(n 5/2 log(n max i,j w ij )) with the algorithm of Orlin and Ahuja [START_REF] Orlin | New scaling algorithms for the assignment and minimim cycle mean problems[END_REF]. Note that, most of the time, algorithms to solve the assignment problem are described as minimum cost assignment with integer costs. However, the maximization cost problem P can be transformed into a minimization cost problem by taking max k, w k -w ij instead of the w ij and costs can be transformed into integers by multiplying them by a suitably large number. In practice, since the obtained problem is a special case of the minimum cost flow problem [START_REF] Ahuja | Network flows: theory, algorithms, and applications[END_REF], it is possible to use any of the algorithms solving it. In particular, it belongs to a class of integer linear problems for which the constraint matrices are unimodular and can be then solved by a linear programming solver (relaxing Constraints 8 from Transforming W to max k, w k -W in order to define a minimization problem, we have:

R ij ∈ {0, 1} to R ij ∈ [0, 1]) knowing that it
W =    
0.3167 0.3167 0.3167 0.3167 0.3167 0.3167 0.2570 0.0559 0.1274 0.3167 0 0.0416 0.3167 0.3167 0.3167 0.3167

    .
The linear assignment problem with cost matrix W admits the following solution:

R =     1 0 0 0 0 0 0 1 0 0 1 0 0 1 0 0     .
Deleting the last line and setting R ij = R ij if w ij ≥ 0 and R ij = 0 otherwise, we get the following final solution:

R =   0 0 0 0 0 0 0 1 0 0 1 0   ,
meaning that objet e 1 is not associated, e 2 is associated with f 4 and e 3 is associated with f 3 .

Finally, we can remark that, while the Hungarian algorithm is historically the most popular method to solve the assignment problem, some refinements and generalizations can be found in the literature (see, for example, [START_REF] Blackman | Design and Analysis of Modern Tracking Systems[END_REF] and [START_REF] Bourgeois | An extension of the Munkres algorithm for the assignment problem to rectangular matrices[END_REF]), in particular to avoid the transformation to square matrices. These algorithms can be used instead of the proposed method.

IV. SENSOR FUSION APPLICATION

In this section, we first show how pairwise mass functions can be computed from object attributes in object recognition applications (Subsection IV-A). We then present experimental results with simulated and real data in Subsections IV-B and IV-C, respectively.

A. Computation of mass functions

Typically, objects are described by a set of attributes. The values of these attributes for each object pair (e i , f j ) can be considered as pieces of evidence regarding the association variable R ij . This evidence may be represented as a mass function m ij in different ways, depending on the nature of the attributes. As an illustration, we will consider three attributes commonly used in object recognition applications: position, velocity and class.

Position: Assume that E and F are sets of objects perceived by two sensors, and each sensor also provides an estimated position for each object. Let d ij denote the distance between the estimated positions of e i and f j , computed using some distance measure (like the Euclidean distance or the Mahalanobis distance if each sensor also returns a covariance matrix). How can a mass function m ij be deduced from d ij ?

Here, it is clear that a single object cannot have two distinct positions and, conversely, two objects cannot occupy exactly the same position. Consequently, a small value of d ij supports the hypothesis R ij = 1, while a large value of d ij supports the hypothesis R ij = 0. Depending on sensor reliability, a fraction of the unit mass should also be assigned to Θ ij = {0, 1}. This line of reasoning justifies a mass function m p ij of the form:

m p ij ({1}) = ρϕ(d ij ) (21a) m p ij ({0}) = ρ (1 -ϕ(d ij )) (21b) m p ij (Θ ij ) = 1 -ρ, (21c)
where ρ ∈ [0, 1] is a degree of confidence in the sensor information and ϕ is a decreasing function taking values in [0, 1]. For instance, the following form may be chosen for ϕ:

ϕ(d) = exp(-γd), (22)
where γ is a positive coefficient. Velocity: Let us now assume that each sensor returns a velocity vector for each object. Let d ij denote the Euclidean distance between the velocity vectors of objects e i and f j . Again, d ij is a piece of evidence about R ij . However, this piece of evidence does not have the same interpretation as the previous one: here, a large value of d ij supports the hypothesis R ij = 0, whereas a small value of d ij does not support specifically R ij = 1 or R ij = 0, as two distinct objects may have similar velocities. Consequently, the following form of the mass function m v ij induce by d ij seems appropriate:

m v ij ({0}) = ρ 1 -ψ(d ij ) (23a) m v ij (Θ ij ) = 1 -ρ 1 -ψ(d ij ) , (23b)
where, as before, ρ ∈ [0, 1] is a degree of confidence in the sensor information and ψ is a decreasing function taking values in [0, 1]. This function can be chosen to have the same form as [START_REF] Ristic | Global cost of assignment in the TBM framework for association of uncertain ID reports[END_REF], possibly with a different coefficient γ . Class: In many applications, objects are categorized in distinct classes such as pedestrian, car, motorcycle, etc. Let Ω be the set of possible classes, and let m i and m j denote mass functions representing evidence about the class membership of objects e i and f j . Such mass functions may be provided, e.g., by evidential pattern classifiers such as described in [START_REF] Denoeux | A k-nearest neighbor classification rule based on Dempster-Shafer theory[END_REF], [START_REF]A neural network classifier based on Dempster-Shafer theory[END_REF]. Ristic and Smets [START_REF] Ristic | Global cost of assignment in the TBM framework for association of uncertain ID reports[END_REF] have considered the object association problem using such class information. However, they made the assumption that equality of class implies equality of objects, a questionable assumption when the number of classes is not much greater than the number of objects. The computation of a mass function m ij on the frame Θ ij from the mass functions m i and m j on Ω can be performed rigorously as follows.

Let S ij denote the event that objects e i and f j belong to the same class, and let Ω ij = {S ij , S ij }, where S ij is the negation of S ij . As shown in [START_REF] Denoeux | EVCLUS: Evidential clustering of proximity data[END_REF], the belief and plausibility of S ij induced by m i and m j have the following expressions:

Bel({S ij }) = ω∈Ω m i ({ω})m j ({ω}) = η ij , ( 24a
) P l({S ij }) = 1 - A∩B=∅ m i (A)m j (B) = 1 -κ ij , ( 24b
)
where κ ij is the degree of conflict ( 6) between m i and m j . The corresponding mass function

µ ij on Ω ij is µ ij ({S ij }) = η ij (25a) µ ij ({S ij }) = κ ij (25b) µ ij (Ω ij ) = 1 -η ij -κ ij . (25c)
Now, it is clear that two objects from different classes cannot be identical, whereas two objects from the same class can be identical or not, which can be formally expressed as follows:

S ij ⇒ (R ij = 0) (26a) S ij ⇒ (R ij = 0) or (R ij = 1). (26b)
Consequently, a mass function m c ij on Θ ij can be computed from µ ij by transferring the mass κ ij to {0} and the remaining mass 1 -κ ij to Θ ij . We thus get:

m c ij ({0}) = κ ij (27a) m c ij (Θ ij ) = 1 -κ ij . ( 27b
)
For each object pair (e i , f j ), a mass function m ij on Θ ij representing all the available evidence about R ij can finally be obtained by combining m p ij , m v ij and m c ij using Dempster's rule:

m ij = m p ij ⊕ m v ij ⊕ m c ij . (28)

B. Simulation experiment

The above approach to object association was tested using simulated data. Each instance of the assignment problem was randomly generated as follows.

a) Experimental settings: We assumed that each of two agents perceives n objects, of which 80% are real objects and 20% are spurious. The position x i of each real object i was generated from a uniform distribution in the square [0, 5] 2 while its velocity v i was generated with direction uniformly distributed in [0, 2π) and with norm uniformly distributed in [0, 0.5]. Objects were assumed to belong to one of two equiprobable classes and to be described by a feature y i normally distributed with standard deviation σ = 2 and means -1 in class one and +1 in class two.

For each real object i, each of the two agents was assumed to get noisy versions of x i , v i and y i defined as follows:

x i = x i + i , v i = v i + i , y i = y i + i , (29)
where i and i are Gaussian noises with mean (0, 0) and variance 0.04I (I being the identity matrix) and i is a Gaussian noise with mean 0 and standard deviation 0.2. Mass functions m p ij and m v ij were computed as described in Subsection IV-A with the Euclidean distance. Each agent was also assumed to compute a mass function m i about the class of object i from its noisy feature y i using the following formula:

m i ({1}) = f 1 ( y i ) f 1 ( y i ) + f 2 ( y i ) , m i ({2}) = 1 -m i ({1}), (30)
where f k is the density distribution of y i in class k.

The spurious objects were generated as real ones, but independently for the two agents. An example of such an association problem with n = 40 objects (including 32 "real" ones) is represented in Figure 1.

b) Performance assessment: The number n of objects was varied between 5 and 80. For each value of n, 30 association problems were generated and solved using the algorithm described in Section III, with ρ = ρ = 0.7, γ = γ = 0.2 and λ = 0. The quality of the association was measured by two criteria: precision and recall defined, respectively, as the fraction of matched pairs that are correct, and the fraction of true object pairs that were matched. Figure 2 shows the mean precision and recall as functions of n in three cases: using the position information only, using position and velocity, and using all three sources of information (position, velocity and class). As expected, the method makes effective use of additional information encoded in the pairwise mass functions, and the performances degrade gracefully with n. Figure 3 shows the mean computing times as a function of n; the association algorithm was programmed in Matlab and run on a Macbook Pro personal computer. c) Influence of ρ and γ: We studied the influence of parameters ρ (= ρ ) and γ (= γ ) on the quality of the association. For each pair (ρ, γ), we computed the F-measure, defined as the harmonic mean between precision and recall:

F = 2 precision • recall precision + recall . (31)
Parameter λ was set to 0. All three sources of information (position, velocity and class) were used. For each pair (ρ, γ), the values of the F-measure were averaged over 30 association problems with n = 20 objects. The results are reported in Table I. As we can see, the choice of parameters ρ and γ does influence the results, and the best performances were obtained for ρ = 0.7 and γ = 0.2. These results show that parameters ρ and γ do have an influence on the quality of the association. If a dataset with matched objects is available, the performances of the system can be significantly enhanced by learning these parameters from the data. We also studied the influence of λ in ( 14) and [START_REF] Kuhn | The Hungarian method for the assignment problem[END_REF]. We recall that this parameter makes it possible to introduce prior knowledge about the number of associations, with positive (respectively, negative) favoring larger (respectively, smaller) numbers of matched object pairs. Consequently, increasing λ can be expected to increase recall (as more true object pairs will be matched) but to decrease precision (as more incorrect pairs will also be matched). An experiment was carried out with n = 40 objects generated as explained above, with λ ranging between -1 and +3. For each value of λ, 200 association problems were generated. Mean running time (in seconds) for our method (maximum plausibility) and Mercier's algorithm, as a function of n.

The results reported in Figure 4 confirm that different values λ indeed result in different trade-offs between precision and recall. However, a higher precision can only be obtained at the cost of a very low recall for that problem. By assigning different weights to precision and recall, it would be possible to find an optimal value of λ. In most applications, however, it might be sufficient to set λ = 0 as a reasonable default value. e) Comparison with Mercier's method: Finally, our algorithm was compared to Mercier's method [START_REF] Mercier | Object association with belief functions, an application with vehicles[END_REF]. However, the comparison could only be carried for small values of n (n ≤ 7) because the computational complexity of this method restricts its applications to very small numbers of objects. On this problem, both methods yield identical solutions in almost all cases. However, as shown in Figure 5, the running time of Mercier's method grows exponentially with n, whereas the time complexity of our method is polynomial.

C. Experiment with real data

Our approach was applied as part of an advanced driver assistance application [START_REF] Fayad | Updating confidence indicators in a multi-sensor pedestrian tracking system[END_REF]. A car was equipped with two sensors: a vision-based Mobileye system and a four-layer Ibeo Alasca-XT Laser scanner. Each sensor has an associated information processing system allowing it to track and classify objects. It is thus a track-to-track association problem [START_REF] Blackman | Design and Analysis of Modern Tracking Systems[END_REF] but we do not use the history of the tracks (or their identity) to match them.

For each sensor, position-based mass functions m p ij were computed using ( 21) and ( 22) with ρ = 0.9, γ = 0.1 and d ij defined as the Mahalanobis distance:

d ij = (x i -x j ) (P i + P j ) -1 (x i -x j ), (32)
where x i denotes the estimated position of the center of object i and P i is the estimated covariance matrix of the estimate. Each sensor also predicts the class of objects, expressed in different frames. The laser classifies the objects as pedestrian or non pedestrian using the algorithm developed by Fayad et al. [START_REF] Fayad | Updating confidence indicators in a multi-sensor pedestrian tracking system[END_REF]; this algorithm computes a mass function over the two classes. The Mobileye system has a finer frame of discernment with five classes: pedestrian, car, truck, motorbike, bicycle. As the latter sensor provides a predicted class, but no confidence in the prediction, a mass 0.9 was assigned to the predicted class, and 0.1 to the whole frame of discernment. For each pair of mass functions m i and m j , a mass function m c ij was computed based on their degree of conflict using [START_REF] Garey | Computers and intractability: a guide to the theory of NP-completeness[END_REF], and it was combined with m p ij using Dempster's rule. We note that velocity information was not used in this application.

The considered dataset contains 58 frames (association problems). The Mobileye sensor detected between 1 and 3 objects (average: 2) while the laser scanner wass less selective and detected between 2 and 23 objects (average: 8.5). An example of an association problem is shown in Figures 6a and6b, with corresponding result shown in Figure 6c. The average precision and recall for our algorithm (with λ = 0) were 0.75 and 0.90, respectively (with standard deviations 0.33 and 0.30). Mercier's algorithm achieved identical results, with considerably longer computing time (119 seconds on average, against 0.23 seconds for our method).

V. THREE-DIMENSIONAL EXTENSION

As an extension of the problem addressed in this paper, we consider in this section the situation where we have more than two sets of objects. To keep the notation simple, and without loss of generality, we will consider the three-dimensional association problem, in which we have three sets of objects. The problem will first be formalized in Subsection V-A and its complexity will then be studied in Subsection V-B.

A. Formalization

Let E = {e 1 , . . . , e n }, F = {f 1 , . . . , f p } and G = {g 1 , . . . , g q } denote the three sets of objects perceived, e.g., by three sensors. We are now searching for three relations R ⊂ E × F , S ⊂ F × G and T ⊂ E × G representing, respectively, the correspondence between objects in sets E and F , F and G, E and G. Each of these three relations has to verify properties (1a)-(1b). Let R, S and T denote the sets of relations R, S and T verifying these properties. Additionally, relations R, S and T cannot be determined independently: for instance, if (e i , f j ) ∈ R and (f j , g k ) ∈ S, then we must have (e i , g k ) ∈ T . More generally, for any three objects (e, f, g) ∈ E × F × G, if any two pairs of objects are related, then the third pair has to be related too. Formally, the three following implications must hold for any triplet (i, j, k):

(e i , f j ) ∈ R and (f j , g k ) ∈ S ⇒ (e i , g k ) ∈ T (33a) (e i , f j ) ∈ R and (e i , g k ) ∈ T ⇒ (f j , g k ) ∈ S (33b) (f j , g k ) ∈ S and (e i , g k ) ∈ T ⇒ (e i , f j ) ∈ R. ( 33c
)
The set of solutions to the association problem is thus the set U of triplets (R, S, T ) ∈ R × S × T verifying (33).

As before, we shall assume that we receive pieces of evidence regarding the association of any pair of objects in E × F , F × G and E × G, and these pieces evidence are encoded as pairwise mass functions m ij• , m •jk and m i•k , for all i, j and k. As before, these mass functions need to be expressed in the common frame of discernment U before being combined by Dempster's rule. Let U ij• denote the set of triplets (R, S, T ) ∈ U such that R ij = 1, and let

m ij• ({1}) = α ij• and m ij• ({0}) = β ij• . To express m ij• in U, we need to transfer the mass α ij• to U ij• , β ij• to U ij• and 1-α ij• -β ij• to U.
The corresponding contour function in U is then defined as follows:

pl ij• (R, S, T ) = 1 -β ij• if R ij = 1, 1 -α ij• otherwise, (34)
for all (R, S, T ) ∈ U, which can be expressed more concisely For each (e i , f j , g k ) ∈ M , we set R ij = S ik = T jk = 1 while all other values defining the relations R, S and T are set to 0. By construction, such relations respect constraints (33) and the value of the objective function is 3r since all the 3r weights associated with the pairs which are kept in the relations have value 1. Now, suppose that there exists a 3-dimensional association (R, S, T ) ∈ U such that the value of the objective function is 3r. Therefore, each of the r objects e i ∈ E is associated with exactly one object f j ∈ F (e i Rf j ) with w ij• = 1 and exactly one objet g k in G (e i T g k ) with w i•k = 1 such that f j Sg k holds and w •jk = 1. For each of the r such associations e i , f j , g k defined by these relations we have (e i , f j , g k ) ∈ M . Indeed, w ij• = 1 implies that ∃y ∈ Y = G such that (e i , f j , y) ∈ M , w i•j = 1 implies that ∃x ∈ X = F such that (e i , x, g k ) ∈ M , w •jk = 1 implies that ∃w ∈ W = E such that (w, e i , g k ) ∈ M , and finally it implies that (e i , f j , g k ) ∈ M since M is pairwise consistent. Thus, there exists a subset M ⊆ M such that |M | = r and no two elements of M agree in any coordinate.

Thus, it is unlikely to solve the 3-dimensional association problem in polynomial time. As future research, we can adapt some methods of the literature for 3DAP to solve our problem. We can then consider branch-and-bound approaches (see, for example, [START_REF] Balas | An algorithm for the three-index assignment problem[END_REF]) or heuristic algorithms (such as described in [START_REF] Pierskalla | The tri-substitution method for the threemultidimensional assignment problem[END_REF]) if the problem is too hard to solve from an operational point of view.

VI. CONCLUSION

Object association is an important problem in a wide range of applications and a key component of many data fusion systems. Belief functions have often been considered as a convenient formalism for representing and combining information in multi-sensor applications (see, e.g., [START_REF] Mercier | Object association with belief functions, an application with vehicles[END_REF], [START_REF] Dallil | Sensor fusion and target tracking using evidential data association[END_REF], [START_REF] Ristic | Global cost of assignment in the TBM framework for association of uncertain ID reports[END_REF], [START_REF] Schubert | Clustering decomposed belief functions using generalized weights of conflict[END_REF]). However, previous attempts to solve the assignment problem within the Dempster-Shafer framework had led until now to ad hoc and sometimes very time-consuming methods (see, e.g., [START_REF] Rombaut | Decision making in data fusion using Dempster-Shafer's theory[END_REF], [START_REF] Gruyer | Credibilistic multisensor fusion for real time application, application to obstacle detection and tracking[END_REF]- [START_REF] Dallil | Sensor fusion and target tracking using evidential data association[END_REF]).

In this paper, evidence about the possible association of any two pairs of objects has been modeled by Dempster-Shafer mass functions defined over the frame of all possible relevant relations between the two object sets. The plausibility of any single relation after pooling all available mass functions can then be computed efficiently and maximized to find the most plausible relation. This problem has been shown to be equivalent to a linear assignment problem, which can be solved in polynomial time using, e.g., the Hungarian algorithm. This solution is thus both optimal and computationally much more efficient than previous approaches to this problem in the belief function framework [START_REF] Mercier | Object association with belief functions, an application with vehicles[END_REF].

As an extension of this work, the three-dimensional association problem has also been considered. This problem occurs, e.g., when fusing reports from three sensors. Although this problem can be formalized in the same way as the two-dimensional one, it was shown to be N P-hard, which makes it unlikely that an optimal solution can be found in polynomial time. The development of heuristic methods to find approximate solutions to this problem is left for further research.

Example 1 .

 1 always has integer solutions. To illustrate the way the object association problem can be transformed into a linear assignment problem, let us consider the following example. Assume that n = 3, p = 4, and the α ij s and β ij s have the following values: the rows and the columns correspond to objects e i and f j , respectively. The corresponding weight matrix W = (w ij ) is:

Fig. 1 .

 1 Fig. 1. Example of an association problem with 40 objects perceived by 2 sensors (a and b). Objects 33 to 40 (marked by a x) are spurious. The size of the circles is proportional to the mass m i ({1}) of class 1.

Fig. 3 .

 3 Fig. 3. Mean running time (in seconds) plus or minus one standard deviation, as a function of n.

Fig. 2 .

 2 Fig. 2. Mean precision (a) and recall (b) as a function of n.

Fig. 4 .

 4 Fig. 4. Mean precision and recall (averaged over 200 association problems with n = 40 objects), as a function of λ.

Fig. 5 .

 5 Fig. 5.Mean running time (in seconds) for our method (maximum plausibility) and Mercier's algorithm, as a function of n.

Fig. 6 .

 6 Fig. 6. Example of an association problem: the laser sensor detects three objects (a) while the Mobileye sensor detects four objects, including a spurious one (b). The association algorithm correctly matches the three real objects (c).

TABLE I F

 I -MEASURE FOR DIFFERENT VALUES OF ρ AND γ (AVERAGES OVER 30 ASSOCIATION PROBLEMS WITH n = 20 OBJECTS, TAKING INTO ACCOUNT POSITION, VELOCITY AND CLASS).

					γ
			0.1	0.2	0.3	0.4	0.5
		0.5	0.641 0.559 0.490 0.425 0.326
		0.6	0.836 0.842 0.744 0.666 0.587
	ρ	0.7	0.795 0.858	0.847 0.797 0.722
		0.8	0.770 0.829 0.846 0.838 0.798
		0.9	0.750 0.804 0.826 0.827 0.821
	d) Influence of λ:

The authors are with Heudiasyc (UMR 7253), Université de Technologie de Compiègne and CNRS, Compiègne, France.

ACKNOWLEDGMENT

This work was carried out in the framework of the Labex MS2T, which was funded by the French Government, through the program "Investments for the future" managed by the National Agency for Research (Reference ANR-11-IDEX-0004-02). The authors would like to thank G. Dherbomez for his help with the acquisition and preprocessing of the experimental data used in Section IV-C, and Vincent Frémont who developed the tools that allowed us to exploit the data from the Mobileye system.

as follows:

Similarly, with obvious notations, mass functions m •jk and m i•k induce the following contour functions in U:

Combining the np+pq +nq pieces of evidence by Dempster's rule yields the following contour function:

Taking the logarithm and, as before, assuming the αs and βs to be strictly smaller than one, we get:

where C is a constant and

The most plausible association (R * , S * , T * ) can thus be found by solving the following binary linear programming problem:

where constraints (41d)-(41f) ensure property (33).

In the special case where all mass functions in one of the three sets {m ij• }, {m •jk } or {m i•k } are vacuous, the problem becomes very simple. For instance, assume that all mass functions m i•k are vacuous. This is the case when, for instance, E, F and G are sets of objects perceived by a single sensor at successive time frames, and we only compute the similarity between objects perceived at two consecutive times. We then have w i•k = 0 for all i and k and the objective function (40) becomes a function of R and S only. As constraints (41d)-(41f) can always be satisfied for some relation T , given R and S, the objective function can be maximized with respect to R and S separately. In that special case, the three-dimensional problem can thus be solved by solving two two-dimensional problems of the form studied in Section III. In the general case, however, the three-dimensional problem is much harder than the two-dimensional one, as will be shown in the next subsection.

B. Complexity analysis

Just as the search for the most plausible relation R * in R is very close to the linear assignment problem, finding the most plausible association (R * , S * , T * ) ∈ U is close to a problem known as the axial 3-dimensional assignment problem (3DAP) [START_REF] Schell | Distribution of a product by several properties, directorate of management analysis[END_REF], which is N P-Hard in the general case. However, several special classes of the problem are polynomially solvable [START_REF] Burkard | Handbook of Combinatorial Optimization -Supplement Volume A[END_REF]. 3DAP can be seen as a special case of our three-dimensional extension with |E| = |F | = |G| = n and where R * , S * ,T * are all bijective relations. It is equivalent to finding exactly n triplets (e, f, g) ∈ E ×F ×G which are pairwise disjoints. The major difference with our problem is that, when searching for the most plausible association (R * , S * , T * ) ∈ U, we can have R * ij = 1 while there is no k such that S * jk = 1 and T * ik = 1. Keep in mind that even if we manage to find a transformation of the most plausible 3-dimensional association problem to 3DAP, it does not mean that we have direct methods to solve our problem (since 3DAP is N P-Hard), neither does it mean that our problem is N P-Hard itself. Unfortunately, it turns out to be the case, as proved in the following proposition.

Proposition 2. The most plausible 3-dimensional association problem is N P-Hard.

Proof: To show this result, we use a reduction of the pairwise consistent 3-dimensional matching problem [START_REF] Garey | Computers and intractability: a guide to the theory of NP-completeness[END_REF] denoted as P C3DM in the following and which is N P-complete: With a given instance of P C3DM , we associate an instance of the 3-dimensional association problem in the following way. We take E = W , F = X and G = Y and then n = p = q = r. For all (e i , f j , g k ) ∈ M , we set w ij• = w i•k = w •jk = 1, while all other weights are set to 0. Note that this reduction is polynomial. We consider the decision version of the 3-dimensional association problem, in which the question is whether there exists a 3-dimensional association (R, S, T ) ∈ U such that the value of the objective function is 3r.

Suppose that there exists a subset M ⊆ M such that |M | = r and no two elements of M agree in any coordinate.


Benoît Lelandais

Su Ruan

Thierry Denoeux

Pierre Vera

Isabelle Gardin

Fusion of multi-tracer PET images for dose painting

Keywords: Dose Painting, Positron Emission Tomography, Information Fusion, Segmentation, Belief Function Theory

Introduction

Medical imaging plays an important role in radiotherapy for treatment planning with the delineation of Gross Tumor Volume (GTV) based on Computed Tomography images. A promising approach to increasing local control of the tumor without increasing normal tissue exposure is to take advantage of intra-tumor heterogeneity using functional imaging. Hall has proposed the concept of dose painting, which involves application of a dose prescription based on biological characteristics of the tumor [START_REF] Hall | Dose-painting by numbers: a feasible approach?[END_REF]. To this end, Positron Emission Tomography (PET) is used for sub-volumes definition, also called Biological Target Volumes (BTV) [START_REF] Ling | Towards multidimensional radiotherapy (MD-CRT): Biological imaging and biological conformality[END_REF]. In the last decade, technological innovations, such as the invention of Intensity Modulated Radiation Therapy (IMRT), made a highly precise radiation application possible [START_REF] Thorwarth | Radiotherapy treatment planning based on functional PET/CT imaging data[END_REF]. As a consequence, high conformal radiation doses can now be given to the tumor while sparing organs at risk and normal tissues. IMRT is based on multileaf collimators leading to a radiation with a spatial resolution in the same order of magnitude than PET images (few millimetres (5 -7 mm)). Dose painting is the subject of a lot of research on its technical feasibility [START_REF] Alber | On biologically conformal boost dose optimization[END_REF], the potential benefit in term of dose escalation within the tumor [START_REF] Thorwarth | Hypoxia dose painting by numbers: a planning study[END_REF], respect for dose constraints to organs at risk [START_REF] Choi | Planning study for available dose of hypoxic tumor volume using fluorine-18-labeled fluoromisonidazole positron emission tomography for treatment of the head and neck cancer[END_REF], as clinical validation [START_REF] Duprez | Adaptive dose painting by numbers for head-and-neck cancer[END_REF].

The main radio-tracer proposed in the literature for dose painting is 18 F-FluoroDesoxyGlucose (FDG), due to the enhanced glucose metabolism of tumor cells corresponding to a high uptake, also called FDG positive tissue. Several other new radio-tracers are of interest, such as 18 F-FLuoroThymidine (FLT) to visualize areas corresponding to a high level of cell proliferation [START_REF] Yang | Imaging of proliferation with 18 F-FLT PET/CT versus 18 F-FDG PET/CT in nonsmall-cell lung cancer[END_REF], and 18 F-FluoroMisonidazol (FMiso) to visualize tumoral regions that have an inadequate supply of oxygen (hypoxic regions) and are radioresistant [START_REF] Chang | A robotic system for 18 F-FMISO PET-guided intratumoral pO 2 measurements[END_REF]. Dose prescription could possibly be adapted in accordance with the radio-tracers. For example, an increase of dose frequency deliverance could be planned for high proliferative cells (i.e., FLT positive tissue), or the dose can be boosted for FMiso positive tissues. Thus, the definition of BTV in each tracer image and their fusion is important to make the radiotherapy more efficient.

Although many methods have been proposed for segmenting positive tissue on PET images in the literature, there is no consensus on this issue.

Most of the methods are proposed to delineate only FDG-PET positive tissues. Example include thresholding methods [START_REF] Erdi | Segmentation of lung lesion volume by adaptive positron emission tomography image thresholding[END_REF][START_REF] Daisne | Evaluation of a multimodality image (CT, MRI and PET) coregistration procedure on phantom and head and neck cancer patients: accuracy, reproducibility and consistency[END_REF][START_REF] Black | Defining a radiotherapy target with positron emission tomography[END_REF][START_REF] Nestle | Comparison of different methods for delineation of 18 F-FDG PET-positive tissue for target volume definition in radiotherapy of patients with non-Small cell lung cancer[END_REF][START_REF] Vauclin | Development of a generic thresholding algorithm for the delineation of 18 FDG-PET-positive tissue: application to the comparison of three thresholding models[END_REF] and region growing methods [START_REF] Haiying | Automatic hot spot detection and segmentation in whole body FDG-PET images[END_REF][START_REF] Krohn | 3D volume and SUV analysis of oncological PET studies: a voxel-based image processing tool with NSCLC as example[END_REF][START_REF] Potesil | Automated tumour delineation using joint pet/ct information[END_REF][START_REF] Green | Semiautomatic volume of interest drawing for (18)F-FDG image analysis-method and preliminary results[END_REF]. The disadvantage of these methods is their sensitivity to noise and heterogeneous uptake of positive tissues. Watershed algorithms are also proposed in [START_REF] Riddell | The watershed algorithm: a method to segment noisy pet transmission images[END_REF][START_REF] Tylski | 18F-FDG PET images segmentation using morphological watershed: a phantom study[END_REF][START_REF] Geets | A gradientbased method for segmenting FDG-PET images: methodology and validation[END_REF]. However, they are sensitive to noise because they are based on gradient measure. Some methods based on probabilistic measures [START_REF] Aristophanous | A Gaussian mixture model for definition of lung tumor volumes in positron emission tomography[END_REF][START_REF] Hatt | A fuzzy locally adaptive bayesian segmentation approach for volume determination in PET[END_REF] and fuzzy measures [START_REF] Zaidi | Fuzzy clustering-based segmented attenuation correction in whole-body PET imaging[END_REF][START_REF] Dewalle-Vignion | A New Method for Volume Segmentation of PET Images, Based on Possibility Theory[END_REF] have recently been proposed. Their advantages are the capability of dealing with noise and/or partial volume effect which is due to the low spatial resolution of acquisition system and the post-filtering applied on PET images. However, in the face of the problems due to heterogeneous uptake of positive tissues, these methods fail. The segmentation of FLT and FMiso-PET images is more complicated than FDG-PET images (see Figure 1). In [START_REF] Hatt | Reproducibility of 18 F-FDG and 3'-deoxy-3'-18 F-fluorothymidine PET tumor volume measurements[END_REF], the authors propose to use their probabilistic method [START_REF] Hatt | A fuzzy locally adaptive bayesian segmentation approach for volume determination in PET[END_REF] for the segmentation of FLT-PET images. For FMiso-PET positive tissue segmentation, a statistical approach based on artificial ants moving in a three dimensional image space has been proposed [START_REF] Haase | Selfreproduction versus transition rules in ant colonies for medical volume segmentation[END_REF]. This specific method has been proposed to address low contrast images. To date, there is no method reported to address more than three tracers. Such a method is challenging to develop due to the partial volume effect, the possibly high noise and the low contrast in PET images. In this paper, we provide a tool to fuse multi-tracer PET images for achieving a good dose painting.

As probability and possibility theories, Belief Function Theory (BFT) is a theoretical framework for reasoning with partial knowledge [START_REF] Dempster | Upper and lower probabilities induced by a multivalued mapping[END_REF][START_REF] Shafer | A mathematical theory of evidence[END_REF][START_REF] Smets | The transferable belief model[END_REF]. Moreover, it has the advantage to contain powerful combination operators [START_REF] Smets | Belief functions: The disjunctive rule of combination and the generalized bayesian theorem[END_REF][START_REF] Denoeux | Conjunctive and disjunctive combination of belief functions induced by nondistinct bodies of evidence[END_REF]. In previous work [START_REF] Lelandais | Segmentation of biological target volumes on multi-tracer pet images based on information fusion for achieving dose painting in radiotherapy[END_REF], we have shown that BFT is a good tool to address multi PET tracers. Our method used BFT to address the imperfections by combining neighboring voxels in a Fuzzy C-Means (FCM) algorithm and to combine multi-tracer PET images. In this paper, we propose a new method that is also based on BFT. The belief masses relative to the tumor are estimated by an extension of the Evidential C-Means (ECM) algorithm which is widely used in data classification to address uncertainty and imprecision. Spatial information is added to the classical ECM to reduce uncertainty due to noise. A new objective function is proposed to design a new algorithm named Spatial Evidential C-Means (SECM) algorithm. Our method is divided into three steps: 1) modeling of uncertainty and imprecision in each mono-tracer PET-image, using spatial information. 2) Reduction of imprecision and uncertainty by combining a measure of partial volume effect specific to the acquisition system and neighborhood information. 3) Fusion of multi-tracer PET images to build parametric maps that describe the locations of glucose metabolism, cell proliferation and hypoxia of tumor.

The manuscript is organized as follow. The proposed framework is presented in the next section. A brief presentation of the BFT is also given in this section. The results on PET images are given in section 3 followed by a conclusion.

Material and methods

As our method is based on belief function theory, we first present it in the following section. Then, the method is presented.

Belief Function Theory 2.1.1. Concept

As probability and possibility theories, belief function theory is a theoretical framework for reasoning with partial and unreliable information. Let Ω = {ω 1 , ω 2 , . . . , ω c } be a finite set of classes, called the frame of discernment. Partial knowledge obtained from information sources is taken into account by assigning Basic Belief Assignments (BBAs) or masses to different subsets of the frame of discernment. Function m is defined as a mapping from

2 Ω = {∅, {ω 1 }, {ω 2 }, . . . , {ω c }, {ω 1 , ω 2 }, . . . , Ω} to [0,1] verifying: A⊆Ω m(A) = 1.
(1)

A hypothesis A is called a focal set if m(A) ≠ 0. The mass m(Ω) represents the degree of ignorance (m(Ω) = 1 corresponds to a total ignorance) and the mass m(∅) is used for representing the conflict.

Using this framework, one can distinguish singletons, such as {ω 1 }, from disjunctions, such as {ω 1 , ω 2 } corresponding to a set of exclusive classes. Singletons allow us to represent uncertainty as in probability theory. A BBA is considered as uncertain when masses are approximately equally distributed over singletons, corresponding to equiprobability in probability theory. Disjunctions allow us to represent imprecision. A BBA is considered as imprecise when non-empty masses are assigned over disjunctions, implying that we do not take part in favor of one class or another.

Basic Belief Assignment (BBA) transformation

From a BBA one can determine a belief function (bel) and a plausibility function (pl) such as:

bel ∶ 2 Ω → [0, 1] A → ∑ B⊆A,B≠∅ m(B), (2)
pl ∶ 2 Ω → [0, 1] A → ∑ B∩A≠∅ m(B).
(3)

These two functions can be used in the last stage of a process of decisionmaking. For decision making, one can also use the pignistic probability function (BetP ) [START_REF] Smets | The transferable belief model[END_REF] defined by:

BetP ∶ 2 Ω → [0, 1] A → ∑ B⊆Ω A∩B B m(B) 1-m(∅) . (4)
where . denotes the cardinality of the considered set. The idea is to convert a BBA into a probability function by redistributing the masses on disjunctions towards singletons.

Combination of BBAs

Two BBAs coming from two information sources and defined on the same frame of discernment Ω can be pooled using combination rules.

Let S 1 and S 2 be two reliable sources of information. The conjunctive combination rule is given by:

m S 1 ∩ S 2 (A) = B∩C=A m 1 (B).m 2 (C).
(5)

It allows us to reduce uncertainty and imprecision by transferring belief masses to conjunctions of hypotheses. When the information contained in sources are in conflict, masses are transferred on ∅ in the open world assumption. In order to satisfy m(∅) = 0 (closed world assumption [START_REF] Smets | The transferable belief model[END_REF]), the Dempster's rule is used [START_REF] Dempster | Upper and lower probabilities induced by a multivalued mapping[END_REF]:

m S 1 ⊕S 2 (A) = 1 1 -κ B∩C=A m 1 (B).m 2 (C). ( 6
)
where κ is the mass on the empty set in the open world assumption. Dempster's rule involves normalizing the conjunctive rule by the conflict between sources.

Both the conjunctive and Dempster's rules are commutative and associative, and can also be generalized to N sources of information.

Evidential C-Means (ECM)

The ECM algorithm has been introduced in [START_REF] Masson | Ecm : An evidential version of the fuzzy c-means algorithm[END_REF]) for dealing with partial knowledge in the context of BFT framework. Let {x 1 , . . . , x i , . . . , x n } be a collection of vectors in R p describing the n objects. Let c be the desired number of classes. Each cluster is represented by a prototype or a center v k ∈ R p associated with each class ω k . Then, for each subset A j ⊆ Ω, A j ≠ ∅, a centroid vj is calculated as the barycenter of the prototypes associated with the classes in A j :

vj = 1 c j c k=1 s kj v k (7)
with:

s kj = 1 if ω k ∈ A j 0 otherwise (8)
where c j = A j denotes the cardinal of A j . The distance d ij between x i and the focal set A j is defined by:

d 2 ij = x i -vj 2 (9)
Let V denotes a matrix of size (c × p) composed of the coordinates of the cluster centers such that V kq is the qth component of the cluster center v j . ECM looks for a partition matrix M = (m ij ) (A j ≠ ∅, A j ⊆ Ω), called a credal partition, of size (n × 2 c ) and for the matrix V by minimizing the following objective function:

J ECM (M, V ) = n i=1 {j A j ≠∅,A j ⊆Ω} c α j m β ij d 2 ij + n i=1 δ 2 m β i∅ ( 10
)
subject to:

{j A j ≠∅,A j ⊆Ω} m ij + m i∅ = 1 ∀i = 1, n. (11)
where m i∅ denote m i (∅), and where δ is a weighting parameter. The empty set is used for outliers. Coefficient α ≥ 0 controls the degree of penalization of the subsets according to their cardinality (α = 2 by default) and coefficient β > 1 controls the fuzziness of the partition (β = 2 by default).

Proposed method: Spatial Evidential C-Means (SECM)

For deriving a credal partition from image data, we propose adding spatial information to the ECM. Some methods [START_REF] Ahmed | A modified fuzzy C-means algorithm for bias field estimation and segmentation of MRI data[END_REF][START_REF] Chen | Robust image segmentation using FCM with spatial constraints based on new kernel-induced distance measure[END_REF] have been proposed to derive a fuzzy partition in image processing by using neighborhood information. Because the performance of ECM has shown superiority in many cases, we propose to extend it by integrating spatial information. Let Φ(i) = {1, . . . , t, . . . , T } the set of the T neighbors of a voxel i, including i. Let t ∈ Φ(i) be a neighbor of i and let b t be a weight function that varies in [0, 1] according to the spatial distance between t and i. b t controls the neighborhood influence and is defined as follow:

b t = exp -d 2 it 2 ⋅ η 2 (12)
where η is a parameter controlling the influence of neighborhood. For our application, η 2 = FWHM 2 8 log 2, with FWHM (Full Width at Half Maximum) corresponding to the spatial resolution of PET images. b t equals 1 when d it = 0, and decreases as d it increases. For the distance d ij in equation 9, we propose the weighted sum of the distances in Φ(i):

∑ t∈Φ(i) b t d tj .
The objective function to be minimized is then given by:

J SECM (M, V ) = n i=1 {j A j ≠∅,A j ⊆Ω} c α j m β ij t∈Φ(i) b t d 2 tj + n i=1 δ 2 m β i∅ (13) subject to {j A j ≠∅,A j ⊆Ω} m ij + m i∅ = 1 ∀i = 1, n, (14)
Note that M and V have been defined before. Using this objective function, the lower ∑ t∈Φ(i) b t d tj is, the higher the mass m ij = m i (A j ) will be, and vice versa.

To model the imperfection from PET medical images, we suggest modifying the objective function slightly. We assume that all the voxels belong either to the background (ω 1 hypothesis) or to high uptake of radioactivity (ω 2 hypothesis). Thus, with the closed world assumption, m(∅) = 0, and the term corresponding to the outlier detection can be removed. The objective function becomes then:

J SECM (M, V ) = n i=1 {j A j ≠∅,A j ⊆Ω} c α j m β ij t∈Φ(i) b t d 2 tj ( 15
)
subject to:

{j A j ≠∅,A j ⊆Ω} m ij = 1 ∀i = 1, n. (16)
To minimize J SECM , an alternate optimization scheme similar to those detailed in [START_REF] Masson | Ecm : An evidential version of the fuzzy c-means algorithm[END_REF]) can be used. First, we consider that V is fixed. We can show that M can be updated using the following equation:

m ij = c -α β-1 j ∑ t∈Φ(i) b 1 (β-1) t d 2 (β-1) tj -1 ∑ A k ≠∅ c -α β-1 k ∑ t∈Φ(i) b 1 (β-1) t d 2 (β-1) tk -1 (17)
Then, we consider that M is fixed. Let B be a matrix of size (c × p) defined by:

B lq = i t∈Φ i b t x tq A j ∋ω l c α-1 j m β ij , l = 1, c, q = 1, p, (18)
and H a matrix of size (c × c) given by:

H lk = i A j ⊇{ω k ,ω l } c α-2 j m β ij t∈Φ i b t , k, l = 1, c (19)
We can show that the minimization of J SECM with respect to V is an unconstrained optimization problem which can be simplified as:

HV = B. (20)
This equation can be solved using a standard linear system solver. The proof is given in Appendix A. By alternatively computing M using ( 17) and solving V using ( 18), ( 19) and ( 20) in an iterative process, the convergence of the algorithm is guaranteed. The process can be applied irrespective of the number of classes. It allows us to consider neighborhood information in deriving the credal partition from image data. For comparison, if ECM allows us to model both uncertainty and imprecision in a one-dimensional feature space (gray level), SECM models the imperfections in a T -dimensional feature space, where T is the number of neighbors, in order to take into account spatial relationships.

Proposed framework

The proposed framework first processes each PET image using SECM to model both uncertainty and imprecision. Then, imprecision and uncertainty are reduced by a combination of different pieces of information. The fusion of the three PET images is finally performed to obtain a parametric map allowing us to help the radiation oncologist for dose painting purpose. The details of each step are given below.

SECM for each PET image

A tissue is said to be positive if it presents a high uptake with one tracer PET image. Two focal sets are considered: the back-ground and positive tissue, respectively identified by hypotheses {H b } and {H p }. Because the belief function framework is used, a third focal set representing imprecision is also considered: {H b , H p }.

After applying SECM on PET images, basic belief masses are spread over hypotheses {H b }, {H p } and {H b , H p }. As a result, we obtain certain, uncertain and imprecise voxels, defined as follows: It results that:

• Certain voxels are those having a high belief mass in favor of {H b } or {H p }.

• Uncertain voxels are those having a medium belief mass in favor of both {H b } and {H p }.

• Imprecise voxels are those having a high belief mass in favor of {H b , H p }.

In PET images, the uncertainty comes from the noise, while the imprecision comes from two phenomena. The first one corresponds to partial volume effect that is present for voxels at the transition between two regions. The second one corresponds to the potential heterogeneity of positive tissues. Heterogeneity is reflected by the presence of inhomogeneity of gray levels in tumor region. Thus, the SECM algorithm allows us to separate these two problems: uncertainty due to noise and imprecision due to both partial volume effect and uptake heterogeneity.

Imprecision reduction for each PET image

As the partial volume effect is caused by the low spatial resolution of the acquisition system, we propose to measure this effect using PET phantom images. A Jaszczak phantom is used, which contains eight spheres whose volumes vary from 1.02 to 98.16 ml, and eight contrasts ranging from 1.7 to 22.9. Because the same acquisition system is used in both cases, PET phantom images and patient images have the same spatial resolution. Our objective is to reduce {H b , H p } representing the imprecision. We define a coefficient γ which reflects the part of imprecision due to the partial volume effect. This coefficient can be measured through the phantom because the sphere volumes are known a priori. The coefficient γ, corresponding to the part of belief from {H b , H p } that should be transferred to {H p }, is learned for each sphere according to its volume and its contrast with the background measured after applying SECM. The measurement is carried out by applying our SECM method on the phantom PET images, and by searching with a dichotomical process, the values allowing us to measure the true sphere volumes. The measurement result obtained is presented Figure 3. The lower both the contrast and the volume are, the higher γ is, except for the smallest volumes.

The use of γ for reducing imprecision is carried out by applying a reinforcement mechanism. First, we propose the following simple mass function having two focal elements for representing the a priori contextual knowledge:

m γ ({H p }) = γ, m γ ({H b , H p }) = 1 -γ. ( 21
)
where H p is the subset to reinforce. Then, the combination of each voxelsource with this BBA is applied using Dempster's rule:

M i = m i ⊕ m γ . ( 22
)
This process allows us to transfer a part of belief to the subset {H p } corresponding to the tumor in PET, and thus reduce imprecision. The reduction of imprecision is applied on the entire image. However, it has an influence for voxels where partial volume effect occurs (where m({H b , H p }) >> 0).

Uncertainty reduction for each PET image

Noise is important in PET images due to low contrast, a short acquisition time period and/or a low activity of the tracer. As proposed in [START_REF] Capelle-Laize | Evidential segmentation scheme of multi-echo mr images for the detection of brain tumors using neighborhood information[END_REF][START_REF] Zhang | Information fusion using evidence theory for segmentation of medical images[END_REF][START_REF] Makni | Introducing spatial neighbourhood in Evidential C-Means for segmentation of multi-source images: application to prostate multi-parametric MRI[END_REF][START_REF] Lelandais | Segmentation of biological target volumes on multi-tracer pet images based on information fusion for achieving dose painting in radiotherapy[END_REF], a solution for reducing noise in the belief function framework is to combine the neighborhood information.

We propose to combine mass functions from neighboring voxels. First, a discounting is done according to the same coefficient b t as before, by transferring a part of belief on the set {H b , H p }:

M t ({H b }) = b t M t (H b ), M t ({H p }) = b t M t (H p ), M t ({H b , H p }) = 1 -b t + b t M t ({H b , H p }) (23)
The discounting operation allows us to reduce the influence of voxels that are physically far away from i before the conjunctive combination. Then, in order to respect the closed world assumption, we propose to use Dempster's rule:

M i (.) = ⊕ t∈Φ(i) M t (.) (24)
The combination of neighbor voxels allows us to remove uncertainty due to noise, by transferring their belief on the singletons. Note that residual imprecision, after the previous step, is also reduced at this step.

Decision for one PET image

For segmenting one mono-tracer PET image, the plausibility function is computed (eq. 3) according to the two classes, and voxels are classified according to their highest plausibility value. Note that this process is not carried out if multi-tracer PET images are considered and fused. Algorithm 1 gives an overview of the proposed method. Choose randomly 2 initial cluster centers = V 0 t ← 0 repeat Compute M t using ( 17), ( 12), ( 9) and ( 7) Compute H t and B t using ( 19), ( 18), and

M t Solve H t V t = B t until V t -V t-1 < ǫ ;

Fusion for multi-tracer PET images

As each PET image provides different information on tumor, the fusion of multi-tracer PET images allows us to carry out the dose painting more efficiently.

Frame of discernment. We consider that the FDG PET image is the reference. It means that the positive tissue volume with FLT and FMiso PET images correspond to sub-volumes of the positive tissue volume with FDG PET image. Considering the fact that diagnosis of cell malignancy is clinically done on FDG PET images, and those using FLT and FMizo are supplementary information, FDG PET images have been considered as the reference and FLT and FMizo images as sub-volumes of those defined on FDG PET images. Thus, five classes are distinguished:

• {N} corresponds to Normal tissues (i.e., background with FDG).

• {M} corresponds to tissues presenting an important glucose Metabolism (i.e., positive tissue with FDG only).

• {P } corresponds to tissues with an important cell P roliferation (i.e., positive tissue with FDG and FLT).

• {H} corresponds to Hypoxic tissues (i.e., positive tissue with FDG and FMiso).

• {F } for F ull uptake corresponds to tissues presenting the three phenomena (i.e., positive tissue with FDG, FLT and FMiso).

The frame of discernment for the three PET images is now: Ω = {N, M, P, H, F }. For each voxel of the three modalities, belief masses are spread over 32 hypotheses: {∅, {N}, {M}, {P }, . . . , {N, M}, . . . , Ω}.

SECM and imperfection reduction. After SECM and imperfection reduction for each PET image, the focal sets are {H b }, {H p }, and {H b , H p }. Because of the frame of discernment used for multi-tracer PET image fusion, and according to the medical interpretation of PET images, the useful hypotheses forming the frame of discernment are given in Table 1. It is assumed that a low uptake in FDG corresponds to a normal tissue {N}, while a high uptake in FDG corresponds to a tumoral tissue that can possibly corresponds to a high or low uptake in FLT or FMiso. Moreover, a low uptake in FLT or in FMiso does not mean that the tissue is normal.

Fusion of multi-tracer PET images for parametric image construction. For the fusion of multi-tracer PET images, we propose fusing each voxel of the three PET images using the conjunctive combination rule (Eq. 5). This combination allows us to discriminate between hypotheses {M}, {P }, {H} and {F } which are very helpful for dose painting.

Decision. Finally, a labeling process is carried out by maximizing the plausibility function according to hypotheses {N}, {M}, {P }, {H} and {F }. This process yields a parametric image presenting the multiple BTVs.

Datasets 2.5.1. PET phantom images

Our method was evaluated on two physical PET phantoms. The first one presents a level of noise similar to that encountered on the FDG and FLT PET images. The second presents a level of noise similar to that encountered on the FMiso PET image.

The FDG and FLT phantom contains nine spheres whose volume varies between 0.43 mL and 97.3 mL. Moreover, five contrasts varying from 2 to 8 are considered. Because the smaller sphere (0.43 mL) was not visible on the two lower contrasts, and the 0.99 mL sphere was not visible on the lower contrast, these three spheres have not been considered.

The FMiso phantom contains six spheres whose volumes vary between 0.54 mL and 26.2 mL. Three contrasts varying from 2.7 to 3.8 were considered. For comparison, the background activity concentration with the FMiso phantom was 17 times lower than that of the FDG and FLT phantom. As a result, the noise was more important in this second phantom. The smallest spheres (0.54 and 1.15 mL) have not been considered because they are difficult to visualize. For the same reason, the 2.63 mL sphere has not been considered with the smallest contrast (2.7).

Acquisition parameters for these two phantoms are given in Appendix B.

Multi-tracer PET images of patients

PET images of patients using FDG, FLT and FMiso tracers were also used for the evaluation of our method. The images were obtained from three patients suffering from lung cancer in the RTEP4 clinical trial [START_REF] Vera | Simultaneous positron emission tomography (PET) assessment of metabolism with 18 F-fluoro-2-deoxy-d-glucose (FDG), proliferation with 18 F-fluorothymidine (FLT), and hypoxia with 18 fluoro-misonidazole (F-miso) before and during radiotherapy in patients with non-small-cell lung cancer (NSCLC): a pilot study[END_REF]. As tumors were paramediastinal, no respiratory motion were taken into account.

Registration of images was made as follow. First, a registration of CT images, acquired for each PET/CT exam, was performed using mutual information [START_REF] Viola | Alignment by maximization of mutual information[END_REF] with FDG PET/CT image as reference. Then, registration parameters were used for PET image registration. A transverse slice is given in Figure 4 for three radio-tracers (patient 1). Each image was scaled to its own maximum Standard Uptake Value (SUV).

Approaches used for comparison

We compared our method on PET phantom images with two other methods recently proposed. The first method is the Adaptive Thresholding (AT) proposed in [START_REF] Vauclin | Development of a generic thresholding algorithm for the delineation of 18 FDG-PET-positive tissue: application to the comparison of three thresholding models[END_REF] for tumor segmentation in FDG PET imaging and inspired from [START_REF] Daisne | Evaluation of a multimodality image (CT, MRI and PET) coregistration procedure on phantom and head and neck cancer patients: accuracy, reproducibility and consistency[END_REF]. The second method is the Fuzzy Locally Adaptive Bayesian (FLAB) method introduced in [START_REF] Hatt | A fuzzy locally adaptive bayesian segmentation approach for volume determination in PET[END_REF], with the following parameters: two classes were used, the size of the neighborhood was 3 × 3 × 3 and Gaussian distributions were assumed. Both methods have been proposed for the segmentation of PET images.

Region of interest definition

As our aim is not to detect tumors, but to estimate their volume, we assume that the lesions have been previously identified by the clinician and that a box encompassing them has manually or automatically been placed. Note that small changes in the placement of the box or its size have no significant impact on the estimation of tumor volume using our method.

Results

Parameter setting

As mentioned in the previous section, four parameters have to be initialized before applying our method. The first two parameters are α and β, which control, respectively the degree of penalization of the subsets of high cardinality and the fuzziness degree. Because only three focal sets are used for segmenting PET images ({H b }, {H p } and {H b , H p }), α acts as a weight for hypothesis {H b , H p }, and β controls the flattening of the fuzzy transition represented by {H b , H p } hypothesis. Because we do not wish to favor {H b , H p }, α has been set to its default value of 2. Coefficient β has been set to 6. This value was chosen empirically using phantom data. Compared to the default value 2, changing the value to 6 has no influence in segmenting positive tissues of large volume, but has an influence on the segmentation of small volumes, which are more subject to partial volume effect. The third parameter η controls the influence of the neighborhood. It is related to the spatial resolution of images. Because F W HM is a measure of the spatial resolution of our PET images, and equals 7 mm, this measure was used for initializing η, as presented in Section 2.3. The last parameter γ is used for reducing imprecision in PET images. As explained in Section 2.4.5, its value is determined using PET phantom images.

On PET phantom images

The results of the sphere volume estimation for each method is presented in Figures 5 and6. Figure 5 shows the results for the FDG and FLT phantom, and Figure 6 for the FMiso phantom. Blue, green and red colors are used, respectively, for SECM, FLAB and AT methods. From left to right, the measured volumes are given from the smallest to the largest sphere. As an indication, the true volumes are represented by a horizontal line for each sphere. The closer the color bar and the line, the better the volume estimation. The measured volumes of each sphere are also given for each contrast. The lighter bar corresponds to the lowest contrasts, and the darker bar corresponds to the highest contrasts.

As we can see on the FDG and FLT phantom PET images (Figure 5), SECM provides the best estimation of sphere volumes. For the largest spheres (> 11 mL), the mean absolute error is about 1.6 ± 1.8 mL for the SECM method, while it is about 6.4 ± 4.7 mL and 5 ± 7.1 mL for FLAB and AT, respectively. For the smallest spheres (< 4 mL), SECM and FLAB yield comparable results. For these two methods, the errors are due to the high partial volume effect that makes the class barycenter estimation difficult. About AT, for small and for large spheres, the measured volumes increase with respect to the contrast. This dependence on contrast variation shows the poor robustness of the method, although its calibration step uses contrast information. The good performance of our method is due to our strategy of first modeling the imperfections and then reducing them according to the volume of positive tissue and its contrast with the background.

On FMiso phantom PET images of low signal to noise ratio (Figure 6), we can see the limit of our SECM method. When a very important noise is present, and when the contrast is very low (about 2.7), SECM gives the worst results. This is due to the very important noise for 2.7 contrast images that induces the presence of regions around the spheres presenting high uptakes. These regions are considered by our method as regions that could potentially correspond to heterogeneity of the positive tissue to segment. That is why an over-estimation of volumes is observed with our method for the lower SNR. However considering 3.3 and 3.8 contrasts only, SECM is the only one method that is robust to sphere size variations. Indeed, the maximum of mean volume error estimations is +3.2 mL for SECM (obtained for the smallest sphere), is -6.5 mL for FLAB (obtained for the largest sphere) and -7.7 mL for AT (obtained for the largest sphere).

On patient data

The results of the belief mass estimation using SECM (see section 2.4.1) are given in Figure 7 for patient 1 according to the three-tracer PET images. The belief mass estimation is represented for each tracer PET image by a map according to the three hypotheses, namely {H b }, {H p } and {H b , H p }. On each map, the lighter (respectively, the darker) the intensity of the color is, the higher (respectively the lower) the mass of the corresponding hypothesis is. A high mass for hypotheses {H b } or {H p } means that the voxel certainly belongs to the corresponding hypothesis. For positive tissues, the high masses are found for voxels for which both themself and their neighborhood present a high uptake. They are mainly found in the middle of positive tissues for the three tracers. A voxel having a high mass assigned to {H b , H p } means that it corresponds to an imprecise region. For the three tracers, high masses are in favor of this hypothesis at the transition between the background and positive tissue. It corresponds to voxels that are subject to partial volume effect. Furthermore, as we can see on the FDG belief mass estimation (see Figure 7(b)), the voxels identified in the white circle have high masses in favor of {H b , H p }. It is an imprecise region whose voxels and their neighbors have a similar uptake. It seems to be due to the presence of heterogeneity in the tumor. Thus, hypothesis {H b , H p } models both partial volume effect and heterogeneity regions for which it is difficult to make a decision. Medium masses for a voxel on both {H b } and {H p } means that the voxel is uncertain. This phenomenon is mainly present for the FMiso PET image due to high noise inside the positive tissue region.

The result of the imprecision and uncertainty reduction steps is shown in Figure 8 for patient 1 according to the three-tracer PET images. As we can see for the three tracers, low masses are assigned to hypothesis {H b , H p }, showing that the imprecision has well been reduced. Moreover, high masses (close to 1) are assigned to {H b } and {H p }, showing that the uncertainty has also been reduced. The important point to observe is that the heterogeneous region in the upper right of FDG PET positive tissue is now considered as certainly belonging to the tumor.

The result of the fusion of multi-tracer PET images is given in Figure 9 for patient 1. Figure 9(d) shows the result of the pignistic probability transformation of masses obtained after fusion. The pignistic probability maps are given according to {M}, {P }, {H} and {F } hypotheses. As in the previous figure, the lighter the intensity of the color is, the more certain is the voxel belonging to the corresponding hypothesis. Figure 9(f) corresponds to the segmented regions obtained by taking the maximum plausibility among the five classes. The colors are defined as follows: the black region corresponds to {N}, the pink region corresponds to {M}, the blue region corresponds to {P }, the green region corresponds to {H} and the yellow region corresponds to {F }. Thus, the fusion allows us to obtain the maps of the specific regions to be irradiated for dose painting task. As a high FDG uptake corresponds to tumor cells, FDG PET images are considered as the reference. Tumor cells are not necessarily associated to a high proliferation, nor to hypoxic cells. It is also possible that low-FDG uptake is associated to high-FLT or FMiso uptake. With our fusion model, this kind of information can be seen on a map corresponding to conflict between data. As presented in Figure 9(e), this map (on the empty set) is also highlighted. The conflict corresponds to hypoxic tissues and/or tissues presenting a high cell proliferation that are normal according to glucose metabolism (i.e., positive tissues according to FLT and FMiso and the background with FDG). It shows positive tissues that potentially correspond to healthy tissues. In our case, conflict indicates the presence of bones that have a naturally high cell proliferation. This conflict is interesting information for dose painting because it could correspond to an organ at risk that would need a low radiation dose.

The fusion results for the two other patients are given in Figures 10 and 11. In Figure 10, a necrosis in the tumor region is present and only a small region at the right side of the tumor presents a high cell proliferation associated with a hypoxia. In Figure 11, we observe that the tumor presents a high cell proliferation, but only few voxels correspond to hypoxia. These results show that our method is efficient even if the size of the region to segment is small. They have been visually evaluated by experts. The use of multiple information conveyed by multi-tracer PET images could have a major impact on the treatment by radiation therapy. The integration of images into a radiation therapy protocol becomes more difficult as the number of images increases. Our results prove that this task can be performed using a fusion process. Our fusion method may be routinely used in the future.

Conclusion

In radiotherapy planning, the conventional FDG tracer (glucose metabolism) is broadly used for tumor depiction in PET imaging. However, in many complex situations, new tracers such as FLT and FMISO can provide complementary information: cellular proliferation and hypoxia, respectively. To date, the segmentation is mainly performed on FGD PET images. In this paper we have proposed a framework based on belief function theory to fuse multi-PET tracers. Because the neighborhood information is integrated in the SECM algorithm, our method allows us to model both uncertainty due to noise and imprecision due to partial volume effect and heterogeneity. A process for reducing the uncertainty and the imprecision is proposed, based on a measurement of partial volume effect from the acquisition system. With multi-tracer PET images, a final fusion step is performed to obtain parameter maps that represent different phenomena observed from multi-tracer PET images. These maps allow us to help the radiation oncologist for the treatment planning because it is possible to adapt the dose prescription according to the different sub-volume defining the tumor.

Our method for volume estimation analysis has been validated using two PET phantom images by comparison with a Bayesian (FLAB) and a thresholding (AT) method. On the first phantom, the results have shown that our method gives the best results for the largest spheres (> 11 mL), and comparable results with FLAB for the smallest spheres (< 4 mL). On the second phantom with low SNR, the results has shown that our segmentation method can achieve a better estimation of the sphere volumes when the contrast is higher than 2.7. The poor result with the 2.7 contrast is due to the presence of regions around the spheres presenting high uptake due to noise. These regions are considered by our method as regions that could potentially correspond to the positive tissue to segment.

The major interest of our method is that it can be applied either for modeling uncertainty and imprecision, for segmenting BTV by reducing imperfection or for segmenting BTV from multi-tracer PET images by fusing them. The results on patient data show the feasibility of using our method for dose painting. Our method is a valuable help to adapt the dose prescription according to the different sub-volume defining the tumor. This information fusion process from multi-tracer PET images had never been proposed before.

The results were also evaluated using phantom data, which are similar as the SNR found on FDG and FLT, and FMiso PET images. Future work will address the evaluation of our method on heterogeneous simulated phantom data. A dose painting strategy is under validation [START_REF] Aerts | FDG for dose painting: a rational choice[END_REF]. The impact of our method for treatment planning will be evaluated in our feature work.

Using (A.3) and (A.4):

λ i β 1 β-1 = ⎛ ⎝ j 1 c α (β-1) j 1 ∑ t∈Φ(i) b 1 (β-1) t d 2 (β-1) tj ⎞ ⎠ -1 (A.5)
Returning in (A.4), one obtains the necessary condition of optimality for M:

m ij = c -α β-1 j ∑ t∈Φ(i) b 1 (β-1) t d 2 (β-1) tj -1 ∑ A k ≠∅ c -α β-1 k ∑ t∈Φ(i) b 1 (β-1) t d 2 (β-1) tk -1 (A.6)
Let us now consider that M is fixed. The minimization of J SECM with respect to V is an unconstrained optimization problem. The partial derivatives of J SECM with respect to the centers are given by:

∂J SECM ∂v l = n i j c α j m β ij t∈Φ i b t ∂d 2 ij ∂v l (A.7) ∂d 2 ij ∂v l = 2(s lj )(v j -x t ) 1 c j (A.8)
∂(e -(v k -vr) 2 ) ∂v l = -2(v kv r )e -(v k -vr) 2 (A.9)

From (A.7), (A.8) and (A.9), we thus have: .11) Setting these derivatives to zero gives l linear equations in v k which can be written as:

∂J SECM ∂v l = 2 n i A j ≠∅ c α-1 j m β ij t∈Φ i b t s lj (v j -x t ) (A.10) ∂J SECM ∂v l = 2 i j c α-1 j m β i,j t∈Φ i b t s l,j 1 c j k s k,j v k -x t (A
i t∈Φ i b t x t j c α-1 j m β i,j s l,j = k v k i j
c α-2 j m β i,j s l,j s k,j t∈Φ i b t (A.12)

Let B be a matrix of size (c × p) defined by:

Figure 1 :

 1 Figure 1: FDG, FLT and FMiso PET images for one patient.

Figure 2 :

 2 Figure 2: Poposed fusion sheme from multi-tracer PET images.

Figure 3 :

 3 Figure 3: The value γ, learned according to the volume and the contrast of spheres.

Algorithm 1 :

 1 Proposed framework algorithm. Data: {x 1 , . . . , x n }: n voxels of one PET image Result: Segmented image according to 2 classes: H b and H p Parameters: η: weighting exponent for neighborhood α: weighting exponent for cardinality (default value 2) β > 1: weighting exponent (default value 2) ǫ: termination threshold (default value 10 3 ) γ: part of imprecision due to partial volume effect SECM:

  Imprecision reduction: Compute m γ using (21) Compute M using (22), M t and m γ Uncertainty reduction: Compute M using (23) and M Compute M using (24) and M Decision: Choose H b or H p for each voxel maximizing (3)

Figure 4 :

 4 Figure 4: FDG, FLT and FMiso PET images for three patients.

  Figure 5: Histograms of measured volumes according to the actual sphere volumes for SECM (a), FLAB (b) and AT (c) methods. The actual volumes are also identified by horizontal lines. The figure is better viewed in color.

Figure 6 :

 6 Figure 6: Histograms of mean measured volumes and their standard deviation according to the actual sphere volumes for SECM (a), FLAB (b) and AT (c) methods. The actual volumes are also identified by horizontal lines. The figure is better viewed in color.

Figure 7 :

 7 Figure 7: Representation of the three transverse slices of the three radio-tracers and the corresponding belief mass estimations according to hypotheses {H b }, {H p } and {H b , H p } using SECM

Figure 8 :

 8 Figure 8: Representation of the three transverse slices of the three radio-tracers and the corresponding belief mass estimations according to {H b }, {H p } and {H b , H p } hypotheses after uncertainty and imprecision reduction steps.

Figure 9 :

 9 Figure 9: Images showing results of the method of information fusion on multi-tracer PET images. (a), (b) and (c) are the initial images using FDG, FLT and FMiso tracers. (d) are the parametric images applying first the BBA estimation of each ROI of each image and then fusing these BBAs. The final result showing how dose painting can be achieved is presented in (f); it corresponds to the segmented regions.

Figure 10 :

 10 Figure 10: Images showing results of the method of information fusion for the second patient.

Figure 11 :

 11 Figure 11: Images showing results of the method of information fusion for the second patient.

Table 1 :

 1 Hypotheses considered after SECM for multi-modal PET images.

	Image	{H b }	{H p }
	FDG	{N}	{M, P, H, F }
	FLT {N, M, H}	{P, F }
	FMiso {N, M, P }	{H, F }

Appendix A. SECM optimization

First, we consider that V is fixed, and we introduce n Lagrange multipliers λ i to solve the minimization problem with respect to M :

By differentiating the Lagrangian with respect to the m ij and λ i and setting the derivatives to zero, we obtain:

We thus have from (A.2):

l = 1, c, q = 1, p, and H a matrix of size (c × c) given by:

V is solution of the following linear system:

PET phantom acquisition parameters

Appendix B.1. FDG and FLT phantom Data have been acquired on a Jaszczak phantom (Data Spectrum Cor-poration™, Hillsborough, NC, USA). It contains nine spheres whose volume varies between 0.43 mL and 97.3 mL (0.43, 0.99, 2.08, 3.78, 11.6, 19.3, 27.9, 58.1 and 97.3 mL). At the beginning of the acquisition, the activity concentration of the spheres was of 41.6 kBq/mL, whereas the activity concentration of the background was 5.4 kBq/mL. For comparison, similar concentrations are encountered on FDG and FLT PET exams (5.0 MBq/kg). Then, activity has successively been added to the background to make acquisitions with 5 different contrasts (7.7, 6.3, 4.9, 3.4 and 2.0).

PET acquisitions and image reconstruction were obtained on a Biograph Sensation 16 Hi-Rez (SIEMENS Medical Solution, Knoxville, TN, USA) following the protocol used in clinical routine. PET emission data were acquired in 3D mode with a 162 mm axial Field Of View (FOV). After the acquisition, the data were corrected for dead-time, random, scatter and attenuation, and reconstructed with an Attenuation Weighted Ordered Subset Expectation Maximization method (AWOSEM: 4 iterations and 8 subsets) in a 168 × 168 matrix (voxel size of 4.06 × 4.06 × 2.0 mm3. A Gaussian post-filtering was applied with a full width at half maximum of 5 mm. The spatial resolution was 6.8 mm in transverse plane at the center of the FOV.

Appendix B.2. FMiso phantom

Data have been acquired on an IEC NEMA phantom (Data Spectrum Corporation™, Hillsborough, NC, USA). It contains six spheres whose volumes vary between 0.54 mL and 26.2 mL (0.54, 1.15, 2.63, 5.53, 11.6 et 26.2 mL). At the beginning of the acquisition, the activity concentration of the spheres was 1.2 kBq, whereas the activity concentration of the background was 0.31 kBq/mL. The objective is to have low concentrations as those encountered on FMiso PET images. For comparison, the background activity concentration is 17 times lower than that of the FDG and FLT phantom. Then, activity has successively been added to the background to make acquisitions with 3 different contrasts (3.8, 3.3 and 2.7).

Except for acquisition time, the same acquisition and reconstruction parameters were used. Data were acquired in List-mode during 21 minutes, and nine bootstrap samples of 7 minutes each were extracted. 7 minutes is also used for FMiso PET image acquisition time. Due to the low activity concentration in the phantom, and thus the low Signal to Noise Ratio (SNR) in the reconstructed images, the use of nine bootstrap samples allows us to better evaluate the methods with respect to noise, by analyzing reproducibility.