EP2859479A1

EP2859479A1 - Procede de determination de volumes fonctionnels pour la determination de biocinetiques

Info

Publication number: EP2859479A1
Application number: EP13728700.9A
Authority: EP
Inventors: Eric Barat; Claude COMTAT; Thomas Dautremer; Thierry Montagu; Régine TREBOSSEN
Original assignee: Commissariat a lEnergie Atomique et aux Energies Alternatives CEA
Current assignee: Commissariat a lEnergie Atomique et aux Energies Alternatives CEA
Priority date: 2012-06-11
Filing date: 2013-06-04
Publication date: 2015-04-15
Also published as: FR2991794B1; WO2013186090A1; FR2991794A1; US20150134296A1

Abstract

L'invention a pour objet un procédé de détermination de volumes fonctionnels pour la recherche de cinétiques représentant l'évolution de la concentration d'un radio-traceur dans une zone de tissu biologique, le procédé étant appliqué à des composantes spatiales (106) et comprend les étapes suivantes appliquées de manière itérative (111) selon un schéma de Monte-Carlo par chaîne de Markov : une étape de génération (108) d'un ensemble (formule A)composé d'un ensemble de cinétiques candidates associées à des valeurs de probabilité d'apparition de ces cinétiques, ces valeurs dépendant de la concentration lambda de radio-traceur; une étape d'étiquetage (109) durant laquelle pour chaque composante spatiale d'indice k (106) les probabilités de sélection des cinétiques de l'ensemble(formule B)sont pondérées (102) en y introduisant une fonction lambda kappa représentative de la concentration de radio-traceur dans cette composante de manière à obtenir un ensemble de valeurs indicatrices (103), une valeur indicatrice D k désignant la cinétique à laquelle est associée la composante spatiale k; une étape de construction de volumes fonctionnels (110), un volume fonctionnel VF j , étant composé de l'ensemble des composantes spatiales qui partagent la même valeur indicatrice D k .

Description

PROCEDE DE DETERMINATION DE VOLUMES FONCTIONNELS POUR LA DETERMINATION DE BIOCINETIQUES

L'invention concerne un procédé de détermination de volumes fonctionnels pour la tomographie par émission de positons et s'applique notamment aux domaines de le tomographie par émission de positions et de la pharmacocinétique.

La tomographie par émission de positons, habituellement désignée par l'acronyme TEP, est une technologie d'imagerie de médecine nucléaire qui utilise des molécules marquées avec un radionucléide émetteur de positons pour imager in vivo les interactions moléculaires de processus biologiques avec un minimum de perturbations. La TEP est la modalité de référence pour l'imagerie moléculaire en raison de sa sensibilité moléculaire de l'ordre de la femto-mole et du caractère quantitatif des mesures qu'elle permet. L'imagerie TEP est utilisée en particulier pour déterminer les volumes fonctionnels désignés dans la suite de la description par l'acronyme VF. Un volume fonctionnel VF correspond à une zone constituée de tissus présentant des propriétés métaboliques identiques.

L'état de la technique en matière d'estimation de volumes fonctionnels procède indirectement en reconstruisant dans une première étape une série temporelle d'images discrétisées puis en estimant les volumes fonctionnels à partir de cette série d'images. Or la diminution de la dose de radio-traceur, appelé aussi radio-pharmaceutique, contribue à dégrader le rapport signal sur bruit de chacune des images ce qui nuit, dans cette approche indirecte, à la détermination des volumes fonctionnels.

L'utilisation en routine clinique de la TEP se fait encore souvent dans un cadre statique, spatial 2D ou 3D, où la variable temps est ignorée, par manque d'outils permettant une analyse des données dynamiques robuste, automatique et validée. Cette utilisation est significativement sous- optimale. En effet, l'imagerie moléculaire TEP est par essence dynamique dans la mesure où la distribution du radio-traceur injecté dans l'organisme évolue dans le temps en fonction des interactions moléculaires entre le traceur et les tissus. L'objectif de l'interprétation de l'examen peut alors se résumer à une caractérisation quantitative et objective des extensions spatiales présentant des comportements pharmacocinétiques temporels similaires, d'où la notion de volumes fonctionnels. Elle est quantitative dans la mesure où elle doit se prononcer sur l'incertitude associée à la concentration moléculaire contenue dans chaque volume fonctionnel en termes de quantité de radio-traceur. Elle est objective dans la mesure où le nombre, le positionnement et l'extension des volumes fonctionnels ne reposent pas sur la supervision d'un expert ayant défini au préalable une ou plusieurs régions d'intérêt.

L'état de la technique dans l'étude dynamique des biocinétiques par imagerie TEP s'appuie sur une segmentation temporelle de l'information afin d'obtenir une série 3D. Cela signifie que l'acquisition est subdivisée en tranches temporelles reconstruites indépendamment les unes des autres. Bien que permettant un suivi dynamique de la distribution du radio-traceur à partir de la série d'images reconstruites, cette approche présente l'inconvénient que seule une partie des observations est utilisée pour reconstruire une image. Cette perte statistique induit fatalement une détérioration de la qualité de la reconstruction d'autant plus que la tranche temporelle est courte.

La plupart du temps, la détermination de l'extension spatiale du volume fonctionnel est définie a priori, à partir de l'imagerie anatomique. Ensuite, on mesure à partir de la série temporelle d'images la concentration d'activité moyenne à l'intérieur du volume fonctionnel, qui correspond en réalité à un volume anatomique, et ce pour chaque tranche temporelle. Des méthodes existantes permettent de déterminer l'extension du volume fonctionnel directement à partir des images dynamiques. Des techniques de séparations de sources sont également utilisées pour déterminer les processus moléculaires distincts à partir des images reconstruites. Cependant, le nombre de cinétiques est fixé dans ces approches. A. Reader et al reprennent ce principe dans leur article intitulé Joint estimation of dynamic PET images and temporal basis functions using fully 4D ML-EM, Physics in Medicine and Biology, 51 (21 ) :5455-5474, 2006, en utilisant un dictionnaire de fonctions temporelles et cherchent la décomposition dans une base fonctionnelle qui ajuste au mieux le jeu de données. À partir d'un nombre limité de fonctions de base temporelles, typiquement une petite dizaine, on peut a priori décrire par leur combinaison un nombre infini de courbes d'activité temporelle, désignées habituellement par l'acronyme TAC venant de l'expression anglo-saxonne « Time Activity Curve ». Ces méthodes de reconstruction 4D permettent d'améliorer le rapport signal sur bruit des images reconstruites. Toutes ces techniques cherchent à procéder à une classification, désignée aussi en utilisant le mot anglais « clustering », de données spatio-temporel des désintégrations, mais procèdent indirectement pour y parvenir puisqu'une étape de reconstructions d'images voxélisées sur tout l'espace est nécessaire. En effet, ces techniques requièrent la discrétisation de toutes les dimensions spatiales et temporelle, entraînant la manipulation et le stockage de gros volumes de données ainsi que des temps de calcul considérables. D'autre part, cette dernière approche ne permet pas de quantifier sans approximation l'erreur d'estimation.

Des travaux sur la reconstruction en tomographie d'émission sans recourir à la discrétisation de l'espace ont été menés dans un cadre bayésien non paramétrique. Cette approche a été étendue à la reconstruction spatiotemporelle comme présenté dans l'article de D. Fall, E. Barat, C. Comtat, T. Dautremer, T. Montagu, et S. Stute intitulé Continuous space-time reconstruction in 4D PET, IEEE Médical Imaging Conférence, 2011 . Cet article est appelé Fall, Barat et al . dans la suite de la description. Cette méthode adresse le problème de la détermination quantitative de l'erreur de reconstruction des volumes fonctionnels et l'objectif de diminution de la dose injectée au patient. Elle propose en outre un « clustering » spatio-temporel dans un espace continu directement à partir des observations. Cependant, dans cette approche, les TAC sont toutes considérées comme des cinétiques normalisées. Cette modélisation, bien qu'elle autorise la détection automatique d'un nombre quelconque de cinétiques, n'est pas optimale dans le sens où la discrimination entre cinétiques n'utilise que leurs formes et ne prend pas en compte la notion de concentration, c'est-à-dire l'amplitude des TAC exprimée en termes de Becquerel par seconde, cette unité étant notée Bq.s^"1. La limite de cette modélisation vient du modèle retenu pour le clustering des composantes spatiales à une cinétique particulière. En effet, celui-ci met en œuvre un modèle hiérarchique où la distribution aléatoire des cinétiques, représentée par /c₀, est elle-même distribuée suivant un processus de Dirichlet d'arbres de Pôlya. Puis, la distribution des composantes spatio-temporelles notée H est choisie comme un processus de Dirichlet dont la mesure centrale fait intervenir K,_Q. Aucune information relative à la concentration n'est présente dans ce modèle stochastique.

Un but de l'invention est notamment de pallier les inconvénients précités.

A cet effet l'invention a pour objets un procédé de détermination de volumes fonctionnels pour la recherche de cinétiques représentant l'évolution de la concentration d'un radio-traceur dans une zone de tissu biologique. Le procédé est appliqué à des composantes spatiales et comprend les étapes suivantes appliquées de manière itérative selon un schéma de Monte-Carlo par chaîne de Markov :

- une étape de génération d'un ensemble A¾ composé d'un ensemble de cinétiques candidates associées à des valeurs de probabilité d'apparition de ces cinétiques, ces valeurs dépendant de la concentration λ de radio-traceur ;

- une étape d'étiquetage durant laquelle pour chaque composante spatiale d'indice k les probabilités de sélection des cinétiques de l'ensemble Λ¾ sont pondérées en y introduisant une fonction \_k représentative de la concentration de radio- traceur dans cette composante de manière à obtenir un ensemble de valeurs indicatrices, une valeur indicatrice D_k désignant la cinétique à laquelle est associée la composante spatiale k ;

- une étape de construction de volumes fonctionnels, un volume fonctionnel VF_j étant composé de l'ensemble des composantes spatiales qui partagent la même valeur indicatrice

Selon un aspect de l'invention, une étape calcule des moyennes des volumes fonctionnels VFj obtenus au cours du processus itératif.

Dans un mode de réalisation, une étape détermine des estimateurs statistiques sur les volumes fonctionnels VFj obtenus au cours du processus itératif de manière à quantifier l'incertitude d'estimation. La variance des volumes fonctionnels VFj obtenus au cours du processus itératif peut alors être déterminée. Un intervalle crédible peut aussi être déterminé.

A titre d'exemple, l'ensemble λ¾ est obtenu par la détermination de deux ensembles de variables scalaires Vj et Ι ,·, ces variables se rapportant à des composantes cinétiques d'indice j, où Vj représente un coefficient de pondération fixe indépendant de la concentration A et v_} représente une valeur de concentration de référence associée à la cinétique d'indice j.

Selon un aspect de l'invention, Vj est déterminé en utilisant l'expression :

Vj

expression dans laquelle :

- pour tout j, m._j désigne le nombre de composantes de la distribution M allouées à la j^ômc composante de H étant une mesure aléatoire distribuée selon un processus de Dinchlet partiellement dépendent, étant une collection indénombrable de mesures de probabilité aléatoires qui admet comme distribution un processus de Dirichlet local ;

- j_Kn désigne le plus grand indice des composantes de tc^ possédant au moins une composante de H allouée ;

- Beta(, ) représente représente la loi de probabilité bêta.

Dans un mode de réalisation, r,- est déterminé en utilisant l'expression :

Γ,Ι^, ϋ, ,, Ζ ^ ^ ιΤ,) avec

U] (Γ/) ex 1 ( ,_r < Tj < b_j' ) x Π (1 - ^ 1 (d (A_fe) Tj) < )} expressions dans lesquelles :

V_j désigne le paramètre caractéristique de la concentration associée à la cinétique j ;

D désigne l'ensemble des variables D_k qui indiquent à quelle cinétique de la composante spatiale k est associée ;

T⁾ _j = {k < _η : D_k = j\ représente l'ensemble des indices k des composantes spatiales associées à la cinétique j ;

Vj = \ k < „ : D_k > j } représente l'ensemble des indices k des composantes spatiales associées à une cinétique d'indice supérieur à j ;

w désigne l'ensemble des pondérations w_k des composantes spatiales _k dans le mélange Dirichlet _H ; Z désigne l'ensemble des paramètres, moyenne et matrice de covariance, z_k associés à la composante spatiale k du mélange Dirichlet M ;

l~ indique ici une génération aléatoire indépendante pour chaque paramètre 1 ;

ex signifie ici «proportionnel à» ;

!() représente la fonction «porte» telle que t(œndition) = i si condition est vraie, et 0 sinon ;

V_] désigne un coefficient associé à chaque cinétique intervenant dans la pondération de ladite cinétique dans ;

d{) représente la distance euclidienne dans R : d(xi , x₂) = χι -r₂ | où

X\ , X'2€ R ;

¾'· est un seuil choisi a priori qui permet de paramétrer A¾ .

Selon un aspect de l'invention, l'étape d'étiquetage est réalisée :

- en générant indépendamment pour tout fc < «_n :

j*

(D*|V, v, Q^*, r, w. Z, C, T) ∑_Pj,_k ôj {-) avec

p_JÂ ex 1 (p_j (X_k ) > v_k ) max (p_j (X_K ) , ς ) J^"Q* (T, )

{i ; C,=k}

et générant indépendamment pour tout k, tel que

j'

avec

P_j,k oc 1 {p₃ (Xk) > !¾) l'iHX (¾· . ς) dans ces expressions,

κ_η désigne le nombre de composantes spatiales auxquelles au moins un événement de coïncidence est associé ;

K* désigne le nombre total de composantes spatiales du mélange Dirichlet H pour une itération donnée ;

v désigne l'ensemble des coefficients V) intervenant dans les pondérations du mélange de cinétiques A¾ ;

v représente l'ensemble des variables auxiliaires _Vk utilisées dans la représentation du processus de Dirichlet de cinétiques ;

Q* désigne l'ensemble des arbres de Polya Q^* caractérisant les cinétiques ;

Γ désigne l'ensemble des paramètres i j caractéristiques de la concentration associée à la cinétique j ;

C désigne l'ensemble des variables de classification c_t qui associent chaque événement de coïncidence i à une composante spatiale du mélange Dirichlet ;

T représente l'ensemble des temps d'occurrences r, des événements de coïncidence ;

p_{j k} désigne la probabilité d'associer la composante k du mélange Dirichlet H à la cinétique j ;

P_j(X) représente la fonction de poids prise en λ associée à la probabilité de la cinétique j de J ¾ ;

6j () représente la fonction mesure de Dirac telle que d, (f) - 1 si j = f et vaut 0 sinon ;

~' indique ici une génération aléatoire indépendante pour chaque variable de classification D_k ;

max{) représente la fonction «maximum de» ; ζ désigne est un paramètre de seuil choisi de façon arbitraire intervenant dans l'implémentation de l'algorithme de « slice sampiing » du processus de Dirichlet local de cinétiques ;

{i : Q - k} désigne l'ensemble des événements de coïncidences i qui sont associés à la composante spatiale fc ;

/Q. représente la cinétique j i.e. la densité de l'arbre de Po!ya j du mélange Dirichlet local de cinétiques ;

J* désigne le nombre total de composantes du mélange Dirichlet local de cinétiques pour une itération donnée.

Selon un autre aspect de l'invention, la concentration Â_k est obtenue par app

avec :

dans laquelle :

E_1¾ désigne l'espérance mathématique par rapport à la loi de probabilité gaussienne de x_k, Λ'(ί ;

N() représente la loi de probabilité gaussienne.

L'invention a aussi pour objet un dispositif de tomographie par émission de positons mettant en œuvre le procédé décrit précédemment.

Le procédé selon l'invention présente l'avantage déterminant d'accroître les distances entre volumes fonctionnels et de ce fait d'augmenter la qualité de séparation.

D'autre part, le procédé selon l'invention réside sur une approche nouvelle permettant d'établir judicieusement la loi a posteriori de la concentration de chaque cinétique, conférant ainsi au procédé la possibilité de fournir une estimation de l'incertitude sur les paramètres asymptotiques de métabolisation d'un tissu. Avantageusement, le procédé selon l'invention permet de diminuer substantiellement fa dose injectée au patient lors d'une imagerie TEP à quatre dimensions tout en fournissant des informations fonctionnelles à l'échelle moléculaire.

Une utilisation répétée de l'imagerie TEP lors de suivi thérapeutique peut alors être envisagée sans préjudice dosimétrique pour le patient. Une utilisation plus développée en pédiatrie et dans la recherche translationnelle peut être également considérée.

En outre, une diminution de la dose de radio-traceur injectée induit automatiquement une diminution de l'exposition aux rayonnements ionisants du personnel médical participant à la préparation et la réalisation de l'examen ainsi que de l'entourage du patient.

D'autres caractéristiques et avantages de l'invention apparaîtront à l'aide de la description qui suit donnée à titre illustratif et non limitatif, faite en regard des dessins annexés parmi lesquels : la figure 1 illustre schématiquement les différentes étapes du procédé de détermination de volumes fonctionnels ;

- la figure 2 illustre la manière dont les poids et variables auxiliaires de la distribution H des composantes spatiotemporelles sont générés ;

la figure 3 illustre schématiquement la manière dont les paramètres spatiaux utilisés par le procédé sont générés ;

la figure 4 illustre schématiquement une manière de déterminer les variables locales et auxiliaires de l'ensemble

la figure 5 illustre la manière dont les cinétiques sont allouées aux composantes de l'ensemble A¾ ;

la figure 6 montre comment les lieux d'émission sont déterminés et comment les composantes de H sont déterminées ;

la figure 7 présente une manière de calculer la concentration \_k de radio-traceur ;

La figure 1 illustre schématiquement les différentes étapes du procédé de détermination de volumes fonctionnels.

Le procédé a pour objectif la détermination de volumes fonctionnels pour la recherche de cinétiques. Il est appliqué à des mesures présentées sous forme d'événements en coïncidence 106 et comprend une série d'étapes appliquées de manière itérative 111 selon un schéma de Monte-Carlo par chaîne de Markov.

Dans une première étape 108 est généré un ensemble t composé d'un ensemble de cinétiques candidates associé à des valeurs de probabilité d'apparition de ces cinétiques, ces valeurs dépendant de la concentration λ de radio-traceur.

Dans une seconde étape dite d'étiquetage 109, pour chaque composante spatiale d'indice k 106 les probabilités de sélection des cinétiques de l'ensemble Λ¾ sont pondérées 102 en y introduisant une fonction À_k. Cette fonction est représentative de la concentration de radio- traceur de cette composante de manière à obtenir un ensemble de valeurs indicatrices 103. Une valeur indicatrice u_k désignant la cinétique à laquelle est associée la composante spatiale k est déterminée. Dans une troisième étape 1 10, les volumes fonctionnels sont construits, un volume fonctionnel VFj étant composé de l'ensemble des composantes spatiales qui partagent la même valeur indicatrice z¾. La méthode de détermination de volumes fonctionnels décrite ci- après permet de pallier l'indépendance du modèle utilisé par Fall, Barat et ai. vis-à-vis de l'information de concentration. Pour ce faire, la concentration de radio-traceur est estimée en tout point de l'espace par marginalisation, c'est- à-dire par intégration sur le temps de la distribution d'activité spatio- temporelle. Cette concentration estimée est alors considérée comme une co- variable aléatoire endogène du problème. Par le terme endogène, il est souligné que cette co-variable n'est pas extérieure car celle-ci est déterminée en commençant la reconstruction. Muni de cette co-variable notée A, la méthode proposée met en œuvre une modélisation dans la distribution des cinétiques dépendante de la concentration de radio-traceur. Ceci revient à substituer à la distribution A;₀ une collection de distributions f .$ =- { h^'£ : A ( A }. Cette collection de distributions admet alors comme distribution un processus de Dirichlet dépendant.

Dans un mode de réalisation, un processus de Dirichlet local IDP d'arbres de Polya est choisi. Une conséquence est que la distribution H des composantes spatio-temporelles n'est plus générée suivant un processus de Dirichlet mais suivant une distribution appelée dans la suite de la description processus de Dirichlet partiellement dépendant et désigné par l'acronyme pdDP venant de l'expression anglo-saxonne « partially dépendent Dirichlet Process ».

Selon un aspect de l'invention, la construction du pdDP est paramétrée par une fonction ψ_χ qui permet de relier la co-variable de concentration A aux paramètres de la distribution marginale spatiale obtenue par intégration sur le temps de la distribution spatio-temporelle. Le problème général de reconstruction d'une distribution de probabilité G( -) dont uniquement des échantillons d'une projection de celle-ci notée F(-) sont observables se formule dans le cadre bayésien non paramétrique pour les problèmes inverses en utilisant les expressions suivantes :

G ~ G

J? (., t) = 7> (.|x) G (dx, t)

K_J. H IF -' F. for i = (1 ) dans lesquelles :

- n est le nombre d'échantillons observés ;

G(.) est une distribution de probabilité aléatoire distribuée suivant ç, définie sur (X _X T._. σ{Χ) ® σ{Τ)) ;

la distribution G{ ) est estimée à partir des observations distribuées selon F( ) : (Υ, τ)' = {(ΥΙ , ΤΊ), . . . . (Υ„, τ_η)} ;

- (-|x) est la distribution de probabilité supposée connue, appelée projecteur, indicée par et définie sur (y, a(y)} ;

la notation _σ(£) correspond à la sigma-algèbre sur i.

L'objet que l'on cherche à reconstruire est la distribution spatiale d'activité dans le champ de vue du tomographe. En TEP, chaque événement en coïncidence est affecté à une ligne virtuelle LOR, acronyme venant de l'expression anglo-saxonne « Line-Of-Response », joignant les centres des deux détecteurs. Une ligne de réponse LOR d'indice l est caractérisée par les paramètres définissant sa position dans l'espace des détecteurs. Une correspondance est ci-après construite entre les paramètres des LOR et l'indice t de telle sorte que la variable y, corresponde à l'indice de la LOR ayant enregistré la ^mc coïncidence observée. Plus précisément, l'observation y_t désigne les coordonnées de la LOR qui joint les deux détecteurs ayant enregistré les photons en coïncidence issus de l'annihilation des positons ayant eu lieu en _Xi.

Dans ce contexte, les termes de l'équation (1 ) correspondent à l'interprétation suivante :

X c R³ correspond à la région de l'espace de reconstruction ;

T c_. f + correspond à la plage temporelle d'acquisition ;

Yi correspond à l'indice de la LOR observée ;

Ti correspond au temps d'occurrence du ^mc événement observé ;

V ( = m x) pour m = 1, . . . , M est la distribution correspondant au projecteur c'est-à-dire à la probabilité d'enregistrer un événement sur la m^ême LOR de détection pour une émission localisée en x, c'est-à-dire le lieu d'annihilation, pour tout x e X,

Cl ) est la distribution spatio-temporelle des émissions correspondant aux événements enregistrés. L'ensemble des observations est tronqué du fait du champ axial limité de l'imageur.

Afin de prendre en compte ces troncatures dans la reconstruction, une LOR d'indice nul est définie, celle-ci correspondant aux émissions non observables de telle sorte que >" = >^» u {<)}. Il est alors possible de définir la probabilité de projection p* (y - m\x) sur ( *. a(y)) et

(yjx) _ P^A (yix)

P (y|x) =

J_y P* (riylx) T^_{= l} V* (j

La densité spatio-temporelle des lieux d'émissions peut alors s'écrire : expression dans laquelle j_y v* (dyjx) désigne la probabilité d'enregistrer une émission localisée en x en tenant compte de la géométrie du système et des propriétés physiques du système de détection.

L'expression (2) définit formellement le modèle proposé qui permet de prendre en compte la concentration spatiale comme co-variable endogène du problème de reconstruction :

= { 1 : λ€ A} ~ ÎDP {y, ¥T_L (A <¾)■ T₀. Φ) (2)

Dans cette expression,

H est une mesure aléatoire distribuée selon un processus de Dirichlet partiellement dépendant pdDP défini comme

ou

H (^■) =∑ ¾ ,¾. (^■)

k=i (3) δ_ξ {·) la mesure de Dirac en ξ.

,V (m,∑) désigne la loi normale gaussienne multi-variée en dimension p ~ 3 paramétrée par sa moyenne m et sa matrice de covariance ∑, de densité : Le modèle génératif défini dans l'expression (2) et la notion de pdDP de l'expression (3) sont explicités ci-après.

L'expression de M dans (3) est un mélange statistique où w_k correspond au poids de la k^émc composante, z_k aux paramètres spatiaux et Q_k à la cinétique, c'est-à-dire la distribution probabilité temporelle, de la fc^êMe composante spatio-temporelle du mélange. Les distributions pour les «¾ et Z_k sont identiques à celles d'un processus de Dirichlet.

On a w = w'i, ~ GE (0) comme indiqué dans l'article de J, Pitman intitulé Some developments of the Blackwell-MacQueen urn scheme, T. F, et al. editor, Statistics Probability and Game Theory, Papers in honor of David Blackwell, volume 30 of Lecture Notes-Monograph Séries, pages 245-267, Institute of Mathematical Statistics, 1996, tel que :

i«i .= î-i ;

pour tout k > 2, u:_k = ¾■ f!*^¹ (i _(¾) où v_L , v₂. . . . sont identiquement distribués suivant une loi Ueta ( l, û), dont la peut être obtenue en utilisant l'expression suivante :

Γ (fi + b)

fl ^ _υ)

Γ (α) Γ («} dans laquelle le symbole r représente la fonction gamma.

Z - Zi , Z₂ , . . . sont identiquement distribuées suivant la loi a priori G₀ qui peut être par exemple choisie comme un modèle Normal-Wishart Inverse {NIW^ ._μ( ,_∑11) tel que z_k - (m_A.,∑_A) avec :

H_k ^{1 L}~ ^" W(∑_k ¹ «o . («ti ¾) ^X).

Où W désigne la distribution de Wishart dont la densité s'exprime en utilisant l'expression (4) dans laquelle n_u est le degré de liberté de la loi Wishart et∑₀ la matrice d'échelle associée.

C'est ici qu'intervient la concentration spatiale comme co-variable endogène pour la distribution des paramètres cinétiques Q_k, Tout d'abord la distribution spatiale s(-) est introduite comme loi marginale de la distribution spatio-temporelle H. L'expression de S( -) est donnée par :

S (.) = / H (; dQ) =

>M(T) _{k l}

dans laquelle M (T) représente l'ensemble des mesures sur T.

On peut remarquer que, par les définitions de m_k et z_k, s(-) est distribué suivant un processus de Dirichlet de paramètre de concentration 0 et de mesure de base a₀ :

S ~ DP (Θ. G₀ )

Q - Q\ . Q>, . . . sont alors indépendamment distribuées selon :

avec x_k = ψ_λ ^ _k, s) où _φχ (■) prend ses valeurs dans A c R.. Par exemple, la fonction ψ_χ peut être choisie comme étant la valeur de l'espérance de la densité spatiale prise sur la composante k :

On peut remarquer que dans le cas où Q_k est indépendant de la marginale S tel que x_k = nous obtenons :

Q_k\Z_k ~ A₀

expressions dans lesquelles F₀ z_k, Q_k ) - G_Q z_k) K°¹ ^ Q_k Cela implique que H ~ DP (Θ, F₀). La distribution de M est un processus de Dirichlet. Dans cette situation la concentration n'intervient pas en tant que covariable endogène. Cependant le modèle couvre alors la possibilité d'avoir recours à une covariable exogène provenant d'une information de type anatomique par exemple.

Il reste à établir la distribution des A'^A. On définit j j = {K : λ G Λ} une collection indénombrable de mesures de probabilités aléatoires qui admet comme distribution un processus de Dirichlet local noté IDF, défini comme suit Soient : v= Vi, v₂, ... sont identiquement distribués tel que Vj ~ lieta(i,;/),

iid

Q* _ Q\.Q?,, , , . sont identiquement distribués tel que Q) ~ W_L (A Q₀) est un arbre de Pôlya à L-niveaux, de paramètres A et de distribution de base Q₀, r— i"₂, ... sont identiquement distribués tel que Γ,· ~ T₀,

Alors,

Ko ( =∑Pi (A) ½;{·)

i=i

avec p (À) = pi (X),p2 (A) , ... défini par l'expression dans laquelle :

_¾(λ) = ¾(λ)Π(ι- Ηλ))

Kj

OÙ

¾(A)= V¾ll(d(A,r_i}< ) et d (λ, λ') représente la distance euclidienne entre A et A'.

On notera que si AC - {K_Q : A - A} ~ 1DP ( , PT (A, <¾) · T₀. ), alors pour tout λ e Λ, κ$\\ ~ DP (//. PT (A, (¾))■ La loi conditionnelle de K sachant A est un processus de Dirichlet d'arbres de Pôlya.

On notera que dans le cas limite où ψ→∞, on a pour tout A e Λ, κ§ = Ko avec ₀ ~ DP (i, PT (A. Q₀)), par conséquent H ~ DP (0. Go x /¾)

II reste à rappeler la définition d'un arbre de Polya pour compléter le modèle. Nous adoptons le formalisme de Lavine présenté dans son article intitulé Some aspects of Polya tree distributions for statistical modelling, Ann. Statist, 20 :1222-1235, 1992, pour la définition de l'arbre de Polya dyadique. Soit E - {(), l } et E^! le /-produit cartésien E χ ^{■ ■ ■} χ E avec E° _ ø. Soit E^* = uf ₀E^L. Soit B₉ = {T} et soit n = {B₀, B B_M, B₀i , . . .} and arbre binaire séparable de partitions de tel que les ensembles H,_{( /+}, sont obtenus par découpe binaire des ensembles B_£1._,._d et j¾ ii_t;..._ti+i définit les ensembles mesurables. L'indice 00 (respectivement 01) est interprété comme 0 (respectivement 1).

On dit que la mesure aléatoire q sur T admet comme distribution un arbre de Polya de paramètres (A, ¾) que l'on notera Q ~ PT (A ¾), s'il existe des réels positifs A - { ,_: : c c E*} et des variables aléatoires yy = {w ·. c c E" } tel que :

la collection W consiste en des variables aléatoires mutuellement indépendantes,

Vt G E*. W_c ~ Bet.a (a_c0. tt_ci).

V I - 1, 2, . . . et Ve f E^!,

avec des facteurs égaux à w où i - if Α· - i.

Considérons la mesure de probabilité g₀. L'arbre de Polya peut être centré autour de Q„ en prenant :

et _e0 = «ei, Ve G E* Dans ce cas, nous notons ₍ = a_; pour ι = 1.2, . . .. et e e E^l. Ce schéma de construction définit l'arbre de Pôlya canonique. Dans ce contexte, si Q ~ PT (A. Qo), alors E [Q) = Q₀.

La construction de la partition peut être poursuivie à l'infini mais comme il n'est pas envisageable de calculer un arbre infini en pratique, les arbres de Pôlya sont souvent simplifiés en spécifiant les paramètres et en construisant la partition jusqu'à un niveau borné L. Si nous notons f_Q (·) la densité de probabilité de Q ~ ΡΤχ (A Qo) alors :

Vt G B₍ f_Q (t) = 2^T'f_Qc (t) · Q (¾, ..._£, ) La justification d'un spécification partielle de l'arbre de Pôlya peut être étayée par le fait que pour les distributions continues, les paramètres augmentent rapidement en fonction des niveaux de l'arbre. La mise à jour de la loi à partir des données influence alors peu la loi a priori car le nombre d'observations dans les niveaux élevés devient faible vis-à-vis du paramètre a priori.

L'objectif étant d'aboutir à un algorithme de reconstruction de volumes fonctionnels, le procédé selon l'invention peut être muni de variables auxiliaires qui vont rendre efficace la procédure itérative de reconstruction et la mise en œuvre d'une méthode de type MCMC, acronyme venant de l'expression anglo-saxonne « Markov Chain Monte-Carlo ».

Nous introduisons c = c_{, c¾ c_m les variables de classification des données aux composantes de H tel que [Zi . Q,) = (Zc, , Q _t) pour tout 1 < n.

De la même manière, soit D = D D₂, . . . tel que Q_k = Q-_Dk pour tout k.

Les variables auxiliaires pour l'échantillonnage par tranches, désigné habituellement par l'expression slice sampling sont décrit ci-après. Afin d'implémenter un échantillonneur de Gibbs efficace, nous suivons une technique dite de slice sampling comme décrite dans l'article de M. Kalii, J. E. Griffin et S. G. Walker intitulé Slice sampling mixture models, Statistics and Computing, 21 (1 ) :93-105, 2011. Cette approche requiert l'introduction de variables supplémentaires. Soit u ~ ui , u-2 , u_n des variables aléatoires uniformes. La densité jointe de (.¾, ¾,«,;), pour une séquence quelconque (¾.) positive peut être déterminée en utilisant l'expression : dans laquelle est fa distribution uniforme sur ]_a, 61 Nous vérifions aisément que par marginalisation sur «, de la formule précédente, nous obtenons la distribution de ^γ, , T, jw. Z, Q). Nous proposons d'adopter une variable de slice (ξ_Α.) dépendante , c'est-à-dire que pour tout k,

dans laquelle c ]o. i] et est indépendant de «_¾. En pratique, des valeurs de correspondant à l'espérance du poids de la première composante ne possédant pas de donnée allouée offrent des propriétés de mélange satisfaisantes (ç % _{n+g^_l+0 - Muni de ce choix de l'expression suivante est obtenue :

( _FÎ T_{i! ¾} Z, Q) - C¹ ∑

dans laquelle les deux sommes sont finies en vertu de #{j : w_} > ε} < oo pour tout > 0.

De façon analogue, soit _v = n. _¾ ... des variables auxiliaires uniformes et ς e ](}, i|. Alors, pour tout k, on a :

(Q_k,

expression dans laquelle \_k = ψ_χ (z_k, s) et ¾^■ ( _k) sont définis supra. Le modèle (2) est ensuite réécrit en lui adjoignant les variables cachées et auxiliaires :

¾ (C,-)

∞

D,, v_k\V, Γ, w, Z ¾' (t;_fc|0, min (ç, _¾ (À,))) P_j (λ_λ) ¾ (¾)

j=i

w ~ GEM (0) Vj ¾ Beta (l, i )

Q; ¾ PT,, (A Qo) (₅) avec

Pj (½) = ¾ (À_fc) Π,^ (l - Vl ( )} Vj (A_fc) _ V₃ 1 (rf (A_{fc ;} ï₃) < φ) et

À partir de ce modèle génératif, il est possible d'exprimer les différentes lois conditionnelles a posteriori qui vont intervenir dans l'échantillonneur de Gibbs.

L'échantillonnage des lois conditionnelles a posteriori par algorithme de Monte-Carlo par chaînes de Markov CMC est décrit ci-après avec l'appui des figures 2 à 7.

À partir de l'expression (5), l'échantillonneur MCMC génère successivement des échantillons provenant des distributions conditionnelles suivantes : Lieux d'émission / allocation aux composantes de H (χ,θΝν,υ,ζ,ςτ,ο,γ,τ)

Paramètres spatiaux : (Z|C, X)

Poids et variables auxiliaires de H : (w,u|C)

Allocation des cinétiques aux composantes de JCQ : (D|V,v,Q^*,I\w,Z,C,T)

Arbres de Pôlya des cinétiques : (Q*|D,C,T)

Pondérations, variables locales et auxiliaires de JCQ : (V,r,_V|D,w,Z,T)

La figure 2 illustre la manière dont les poids et variables auxiliaires de la distribution H des composantes spatio-temporelles sont générés, c'est- à-dire la manière dont des échantillons de distribution conditionnelle ( .ujC) sont générés.

Nous notons le nombre latent de composantes du mélange spatio-temporel pour n observations et n_À. = #{i _: < = k} pour k < κ_η. Nous adoptons l'approche présenté dans l'article de Kalli et al. Précédemment cité et procédons à la génération conjointe de (w.u|C) en échantillonnant tout d'abord w w₂... ,w_K C, puis «(ωχ, «/_Kr_,c, et finalement w_Kn+2, , .. Wi, w₂, ... , w_Kn,C

Dans ce cadre, en suivant la définition neuf et le corollaire vingt de l'article de J. Pitman intitulé Some developments ofthe Blackwell-MacQueen urn scheme dans T. F. et al., editor, Statistics, Probability and Game Theory ; Papers in honor of David Blackwell, volume 30 of Lecture Notes-Monograph Séries, pages 245-267, Institute of Mathematical Statistics, 1996, nous sommes en mesure d'exprimer la loi a posteriori du processus de Dirichlet comme la somme d'une somme finie de composantes dont les poids sont distribués par une distribution de Dirichlet distribution et d'un processus de Dirichlet indépendant des données. La présence de la variable auxiliaire u permet de ne générer qu'un nombre fini κ^* de poids et de lieux d'émission lors de la génération du processus de Dirichlet indépendant.

Dans un premier temps 200, les w_k sont générés pour t. < κ_η. (w₁,w₂,...,w_Kn,r_Kn\C) ~ Dirichlet (m, n₂ n_Kn1 Θ) Dans un second temps 201 , les échantillons UÎ , W₂, ... c sont générés. ΐΗ |C '~ ^' U {li_t |0, min (w_r, , ζ))

XL* = iiiin(u)

Dans un troisième temps 202, les échantillons «¾. sont générés pour k > «„. Tant que _rfc__J > «*,

U_k ~ Beta (^'! , Θ) 7¾ = r_fc i ( i - ¾)

Dans un quatrième temps 203, il est posé _K" = mm {{k ·. r_k < u*} ).

Clairement, _M¾ < u * pour tout k > κ\ Pour cette raison, un jeu fini de «_¾ est suffisant pour l'allocation des données aux composantes. Ajoutons que nous abandonnons dans les calculs également tous les w_k tel que w_k < υ'- pour k < κ_Η.

La figure 3 illustre schématiquement la manière dont les paramètres spatiaux sont générés, c'est-à-dire la manière dont des échantillons de distribution conditionnelle z C_; X) sont générés.

Dans un premier temps 300 sont générés (Z C. X) pour tout k < n_n, à partir d'une loi Normale-Wishart Inverse dont les paramètres sont mis-à-jour à partir des observations de par la propriété de conjugaison :

avec

Qk = Q + ^""k

ρ μυ + nk fîk

ΰ + n_k

1l_k ∑ ¾

{i:C',=Jt}

et avec

=

Dans un deuxième temps 301, les échantillons Z_fesont générés :

La figure 4 illustre schématiquement une manière de déterminer les variables locales et auxiliaires de l'ensemble A¾. Pour cela, (V, Γ, v|D. w, z , T) doivent être générés.

Nous considérons comme loi a priori r_f, (.) Pour tout j, nous définissons

V_j = {k< κ_η : D_k = j}

V† = (k < κ_η : D_k > ./}

Et désignons pour tout j, ^mj = comme le nombre de composantes de M allouées à la f^mc composante de f $.

Nous définissons également le plus grand indice des composantes de tc possédant au moins une composante de H allouée :

max({j : rn, > U})

Nous adoptons le même schéma de génération que pour l'échantillonnage de (w.u'C), c'est-à-dire que nous générons en premier lieu ( ,r,-|D, ,z,T) pour j<j_Kn 400, puis (v\O,v ... , V_JriH,Ti, , . * ,T_JKU. ,Z) 401 et, finalement, Vj et V_j pour j > j_Kn 402.

Pour tout j < j_Kn, sont générés indépendamment : Pour tout j < J_KII, nous posons aj p -^* inax ( max ( , - Ψ) t «r

puis générons :

OU

U (Γ, ) « 1 ( ,_Γ < Γ, < ¾· _Γ) Χ Π (1 - ¾ 1 (d (X_k, Tj) < T»))

L'échantillonnage de ¾J requiert une étape Metropolis-Hastings (MH). Nous proposons un échantillonneur MH indépendant avec comme distribution candidate :

Ainsi, à l'itération t ι l, pour tout j < j^,

- Générer r* ~ u (r '_;-_,r, ¾_.r)

- Calculer le ratio d'accep

où est la valeur de la chaîne de Markov à l'itération t.

- Tester l'acceptation de la proposition 1^"'* i.e. affecter r†^H) = Γ* avec une probabilité

sinon affecter rf^{+l )} = r . - Pour tout k < K_ilt poser _Pi%k = _ . (1 - Vj i [d (X_k, Γ,·) < f>)) et

| = ιηΐιι(¾,ι,...,^·,_Κι1).

- Pour tout k < «„, générer indépendamment

¾|D, ½,..., j_Kn , Γ_ΐ5... , Tj_Kn , w, Z ¾> W (¾[0, min (p_¾ (X_k) , ç)) et poser

u* = initi(v)

- Générer v₃- et i_j pour j > j_Kn. Tant que p*_, > ¾,^»*,

¾· ~ Beta

- Pour tout k < K,,, poser _Pjsk = _Pj__l>k (l - V₃1 (ά{Χ^Γ₃) < V)J et

} = min(p_j,i,...,^._K.,).

- Poser J* - min ({j : p* <

La figure 5 illustre la manière dont les cinétiques sont allouées aux composantes de x¾, c'est-à-dire comment est généré.

Dans un premier temps 500 est généré indépendamment pour tout k < κ_η :

r

(D_k\ V, v, QM\ w. Z,C.T) ¾?∑ p _k Sj (·) avec

pjj_t oc 1 (pj (X_k) > 1·/·) ma (¾^■ (X_k) . ç) JJ /g; (7*)

{t:C,-ft}

Dans un deuxième temps 501 est généré indépendamment pour tout k, te! que H_n < k < K* : avec

P_j,k « 1 ÎP_j ί½) > ¾) max (pj (A_fc) , ς)

11 faut ensuite générer les arbres de Polya des cinétiques générés.

D'après la propriété de conjugaison de l'arbre de Polya, pour tout j < j_{H n} , on génère indépendamment (Q* D. C, T) suivant une distribution d'arbre de Polya avec des paramètres mis-à-jour A_j' - {<y_J' _t ■ ? e tel que

avec

et

n = # {i€\ l n} : (Dr, = i) A (ï| e fi_t)|

Puis, pour tout j _n < j < j*. générer indépendamment ^ ρ, Ο, Τ ^ Ρΐ (Â. Qo)

La figure 6 montre comment les lieux d'émission sont déterminés et comment les composantes de H sont déterminées. Cela se fait en générant (C. Xjw, u. Z, Q^*. D, Y, T).

La loi conditionnelle présente une forme qui ne permet pas une génération directe. En effet, elle fait appel au projecteur v {Y_l\X_l) qui intègre des caractéristiques géométriques et physiques du dispositif de détection, ce qui exclut en pratique de trouver une expression analytique pour la distribution offrant un espoir de génération simple. Nous avons donc recours à une étape de Metropolis-Hastings. Nous commençons par exprimer la distribution conditionnelle en utilisant la règle de Bayes. , Xi lYi, T u w, Z, Q*, D V x (¾ (Q, ¾) avec fc. (Ci-, Xi) oc m) mox (w_k, Q f_Ar (¾[¾) _Q¾ (T.) ¾ (£¾) et ^(y'iiA^'i) correspond à la vraisemblance du projecteur pour le jeux d'événements observées (tronqué). Comme cette vraisemblance n'admet pas de forme standard permettant l'échantillonnage direct, nous utilisons un algorithme de Metropolis-Hastings indépendant et devons définir une loi candidate que nous notons G^* pour / = i, .. , , n. À chaque itération t + i nous devons générer (C*.X-) ~ G; et calculer le ratio d'acceptation. Cette étape est la plus coûteuse en termes de calculs.

Nous proposons ainsi une loi candidate de la forme :

fct (Ct,XÏ) * f (λΪ ,,,ΣΪ-) x f_ôi iC X )

Nous remplaçons le projecteur ρ(ΐίΐλ',) par M dans la loi conditionnelle pour construire la distribution candidate 600, où m'_y et ∑ sont choisis de façon à approximer convenablement τ(γ^χ_{). Nous pouvons alors exprimer /£(·) en utilisant l'expression :

«*

Je; (C*, Xt) oc∑l(«¾ > _Ui) max (w_k, Q r _xJ (Xt] _k._Y: , _i) (C?)

k=l ^k

avec pour tout k < K* et m < Λί ,

= ((∑, )-' +Σ ¹)^"1

Il est alors possible de procéder à l'échantillonnage de (C,,¾), pour tout i < n à l'itération t + i.

C* est généré indépendamment 601 tel que :

K*

Ci ~∑*'*.t<M-)

fc=l avec

w_k,i oc 1 (w_k > ιι,.) max(¾., )¾F, Q._¾ (Γ,) et E^_!^-L

Ensuite x* est généré indédendemment 802 tel que : X^* '~ A^r (X^* 1 rnr , y, ,∑(? , _f )

Les ratios d'acceptation > uite calculés 603 : _u (< ,*?), <cf> x«>) = ^, _ω (x_t,x )

On pourra remarquer que <¾(¾,„¥,) et c disparaissent du ratio d'acceptation.

11 reste alors à tester 604 l'acceptation de (C*, A^',^*), c'est-à-dire d'affecter (c⁺¹⁾.Af⁺ⁱ⁾) = (C7*, *) avec la probabilité

sinon affecter (c|'^+,),Af^{+ 5}) = (c .A\^(iî}) La figure 7 présente une manière de calculer la concentration \_k.

L'évaluation de x_k est nécessaire pour la génération de ;v.r.v|D.w, Z,T) et de (Ό|ν,ν,<¾',Γ,_ν,ζ,ο,τ>

L'expression de \_k, pour le choix particulier de é_\(Z_k.S) donné précédemment, peut être déterminé d'après l'équation (5) en utilisant l'expression :

Cette expression peut être évaluée par une technique d'échantillonnage d'importance. Pour cela, pour un nombre de particules donné τ_λ et pour tout k, on génère indépendamment 700 pour tout ι < Ί\ : Ensuite, on calcule 701 pour tout ι < Ί :

Dans un troisième temps 702, est évalué :

Reste ensuite à déterminer les volumes fonctionnels VF 105,

fois que l'algorithme MCMC a atteint l'équilibre, l'échantillonneur fournit des tirages de la loi jointe a posteriori de (X, C. D, w, u. z, Q* Y. T). Pour chaque tirage de paramètres retenu, il est possible de construire un échantillon de la distribution d'activité spatiotemporelle des événements enregistrés à l'itération {t) : f f (x. <) = -v (^χΙ¾) ÎQ_k (t)

(6)

La densité spatio-temporelle correspondant à tous les événements, qui ne sont pas nécessairement observables, est obtenue par normalisation.

De la même façon, l'entière distribution des volumes fonctionnels est accessible. Nous marginalisons sur le temps et définissons la distribution spatiale d'un volume fonctionnel j associé à la cinétique Q^* à l'itération (t) en utilisant l'expression :

(7) Tous les éléments sont alors disponibles pour déterminer les estimateurs de distribution d'activité spatio-temporelles et des volumes fonctionnels VF pour T itérations de l'algorithme retenues. Les expressions suivantes peuvent alors être utilisées pour cela : (8) et

*=ι ^' (9)

De la même façon, toute autre fonctionnelle des distributions a posteriori peuvent être calculées, comme une variance conditionnelle ou un intervalle de crédibilité.

Claims

REVENDICATIONS

Procédé de détermination de volumes fonctionnels pour la recherche de cinétiques représentant l'évolution de la concentration d'un radio- traceur dans une zone de tissu biologique, le procédé étant appliqué à des composantes spatiales (106) et comprend les étapes suivantes appliquées de manière itérative (111 ) selon un schéma de Monte- Carlo par chaîne de Markov :

- une étape de génération (108) d'un ensemble λ¾ composé d'un ensemble de cinétiques candidates associées à des valeurs de probabilité d'apparition de ces cinétiques, ces valeurs dépendant de la concentration A de radio-traceur ;

- une étape d'étiquetage (109) durant laquelle pour chaque composante spatiale d'indice k (106) les probabilités de sélection des cinétiques de l'ensemble A¾ sont pondérées (102) en y introduisant une fonction x_k représentative de la concentration de radio-traceur dans cette composante de manière à obtenir un ensemble de valeurs indicatrices (103), une valeur indicatrice D_k désignant la cinétique à laquelle est associée la composante spatiale k ;

- une étape de construction de volumes fonctionnels (110), un volume fonctionnel V F_} étant composé de l'ensemble des composantes spatiales qui partagent la même valeur indicatrice D_k.

Procédé selon la revendication 1 comprenant une étape de calcul des moyennes des volumes fonctionnels VF_j obtenus au cours du processus itératif. 3- Procédé selon l'une des revendications précédentes comprenant une étape de détermination d'estimateurs statistiques sur les volumes fonctionnels V ) obtenus au cours du processus itératif de manière à quantifier l'incertitude d'estimation.

4- Procédé selon la revendication 3 dans lequel est déterminé fa variance des volumes fonctionnels v t) obtenus au cours du processus itératif.

5- Procédé selon l'une des revendications 3 ou 4 dans lequel est déterminé un intervalle crédible.

6- Procédé selon l'une des revendications précédentes dans lequel l'ensemble x¾ est obtenu par la détermination (100) de deux ensembles de variables scalaires v_} et Γ,, ces variables se rapportant à des composantes cinétiques d'indice j, où Vj représente un coefficient de pondération fixe indépendant de la concentration λ et v_} représente une valeur de concentration de référence associée à la cinétique d'indice j.

7- Procédé selon la revendication 6 dans lequel V_j est déterminé en utilisant l'expression :

expression dans laquelle :

- pour tout j, ni_j désigne le nombre de composantes de la distribution H allouées à la f^mc composante de /¾\ H étant une mesure aléatoire distribuée selon un processus de Dirichlet partiellement dépendent, |j étant une collection indénombrable de mesures de probabilité aléatoires qui admet comme distribution un processus de Dirichlet focal ;

- j_lin désigne le plus grand indice des composantes de /c possédant au moins une composante de M allouée ;

- Beta(, ) représente représente la loi de probabilité bêta.

- Procédé selon l'une des revendications 6 ou 7 dans lequel i , est déterminé en utilisant l'expression :

avec

u* (Γ,) « i ( _Γ < Yj < ¾._r) x H (1 - v_d i (d (X_k, Tj) < φ)) expressions dans lesquelles :

T_j désigne le paramètre caractéristique de la concentration associée à la cinétique j ;

D désigne l'ensemble des variables u_k qui indiquent à quelle cinétique de ¾ la composante spatiale k est associée ;

Vj - {A- < Kn■ D_k - j} représente l'ensemble des indices k des composantes spatiales associées à la cinétique j ;

V = {k < κ„ ; l^' > j} représente l'ensemble des indices k des composantes spatiales associées à une cinétique dindice supérieur à j ;

w désigne l'ensemble des pondérations w_k des composantes spatiales _k dans le mélange Dirichlet _H ; z désigne l'ensemble des paramètres, moyenne et matrice de covariance, z_k associés à la composante spatiale k du mélange Dirichlet H ; indique ici une génération aléatoire indépendante pour chaque paramètre Γ, ;

ex signifie ici «proportionnel à» ;

l() représente la fonction «porte» telle que i(eonditt<m ) - 1 si condition est vraie, et 0 sinon ;

Vj désigne un coefficient associé à chaque cinétique intervenant dans la pondération de ladite cinétique dans Λ¾ ;

d{ ) représente la distance euclidienne dans R : d(x₁ , x₂) = ι - x₂ | où x\ . x- € R !

^■ est un seuil choisi a priori qui permet de paramétrer .

Procédé selon l'une des revendications précédentes dans lequel l'étape d'étiquetage (109) est réalisée : en générant (500) indépendamment pour tout k < n_n : avec

V₃.K oc 1 (p_j (A*) > i _k ) max ( j (X_k ) . ç) JJ /Q* (T. ) et

en générant (501 ) indépendamment pour tout k, tel que „ < k < * : avec

Pj,k oc t Û¾ (A_fc) > ·¾) max (j¾ (A_FC) , c) et ∑ :_{: ¾fc} = i- dans ces expressions,

V désigne l'ensemble des coefficients Y_j intervenant dans les pondérations du mélange de cinétiques A¾ ;

q* désigne l'ensemble des arbres de Polya Q* caractérisant les cinétiques ;

Γ désigne l'ensemble des paramètres ij caractéristiques de la concentration associée à la cinétique j ;

c désigne l'ensemble des variables de classification c_t qui associent chaque événement de coïncidence ί à une composante spatiale du mélange Dirichlet H ;

T représente l'ensemble des temps d'occurrences r_t des événements de coïncidence ;

p_{j k} désigne la probabilité d'associer la composante k du mélange Dirichlet H à fa cinétique j ;

P_j iX) représente la fonction de poids prise en λ associée à la probabilité de la cinétique j de j¾ ;

6_j() représente la fonction mesure de Dirac telle que <¾{/) - 1 si j = f et vaut 0 sinon ;

'~ indique ici une génération aléatoire indépendante pour chaque variable de classification n_k ;

ma.\() représente la fonction «maximum de» ; ( désigne est un paramètre de seuil choisi de façon arbitraire intervenant dans l'implémentation de l'algorithme de sliœ sampling du processus de Dirichlet local de cinétiques ;

{i : C: = k} désigne l'ensemble des événements de coïncidences i qui sont associés à la composante spatiale k ;

f_Q- représente la cinétique j Le. la densité de l'arbre de Po!ya j du mélange Dirichlet local de cinétiques ;

10- Procédé selon l'une des revendication précédentes dans laquelle la concentration À_k est obtenue ar a lication d'une fonction φ_λ définie telle que = Φχ

dans laquelle :

E_îk désigne l'espérance mathématique par rapport à la loi de probabilité gaussienne de x_k, Ai(.r_k \z_k) ;

N() représente la loi de probabilité gaussienne.

1 1 - Dispositif de tomographie par émission de positons mettant en œuvre le procédé selon l'une des revendications 1 à 10.