Brevet FR3051937A1 — publié le 20171201

PROCEDE ET DISPOSITIF DE CARACTERISATION D'UN FAISCEAU DE CORPUSCULES L'objet de l'invention concerne un procédé et un dispositif permettant la caractérisation d'un faisceau de particules, au sein d'un accélérateur de particules, par exemple, afin notamment de déterminer la dose déposée dans un milieu par les corpuscules émis. Elle s'applique, par exemple pour caractériser les paramètres physiques d'un faisceau d'irradiation de rayons X ou d'électrons modélisés par une simulation de type Monte Carlo. Elle est utilisée, par exemple, dans le domaine de la radiothérapie, de l'imagerie médicale utilisant les rayons X afin de calculer les doses utilisées. Les logiciels de calcul de la dose reçue par un patient au cours d'un traitement de radiothérapie ou lors d'un examen d'imagerie par rayons X nécessitent une connaissance des propriétés physiques du faisceau utilisé pour l'irradiation. Par propriétés physiques, on entend l'énergie, la direction et la position des particules présentes dans ce faisceau. Le calcul de la dose délivrée repose de plus en plus sur des méthodes de calcul de type Monte-Carlo, MC, en raison de leur supériorité à délivrer une information précise, quelques que soient les situations cliniques rencontrées. Avec ces méthodes MC, on réalise d'ordinaire le calcul de dose en deux étapes : la première consiste à calculer un fichier d'espace des phases (noté PSF), qui contient la description physique du faisceau d'irradiation délivré par le dispositif médical thérapeutique, en radiothérapie, ou diagnostique en imagerie médicale par rayons X. Dans une seconde étape, on utilise cet espace des phases pour réaliser le calcul de dose proprement dit, dans une scène géométrique qui décrit par exemple le patient à soigner. Dans une approche Monte Carlo, le faisceau d'irradiation est représenté par un échantillon aléatoire, i.e., une collection de particules générées suivant la distribution de probabilité caractérisant le faisceau. Le problème technique posé porte sur le fichier d'espace des phases. Ce fichier 3 0 5 1 9 3 7 2 stocke toutes les particules appartenant au faisceau et possède la structure suivante : les propriétés de chaque particule du faisceau sont stockées sur une ligne. Ce choix de représentation entraîne trois limitations techniques : 1 - la taille du fichier : pour que le calcul de dose par méthode MC soit 5 représentatif sur un plan statistique, il est nécessaire de stocker un nombre d'échantillons important, par exemple plusieurs millions. Le fichier peut donc atteindre plusieurs centaines de giga-octets, voire plusieurs téraoctets, et devient donc difficilement gérable lorsqu'il est intégré à un logiciel de calcul de dose, 10 2 - le nombre limité de particules stockées dans le fichier peut être à l'origine d'erreurs systématiques, à l'image d'un bruit, sur la représentation du faisceau. Ces erreurs se propagent dans le calcul de dose ultérieur ou de biais si le fichier est lu et utilisé trop souvent lors du calcul de dose, 3 - le caractère non flexible de ce mode de représentation : le fichier 15 d'espace des phases décrit l'état du faisceau à un instant donné et correspond à un paramétrage donné. Si le paramétrage du faisceau change, même légèrement, le fichier doit être de nouveau calculé. Ce mode de représentation est donc peu propice pour prendre en compte les légères variations physiques existantes entre deux faisceaux issus du même 20 dispositif d'irradiation. Il existe donc un besoin de fournir un procédé et un système permettant de stocker de façon plus compacte les caractéristiques physiques d'un faisceau d'irradiation, de générer le nombre de particules nécessaires au calcul de dose, sans limite intrinsèque, de manière flexible et rapide. 25 L'article de C.M. Ma, B.A. Faddegon, D.W.O. Rogers et T.R. Mackie intitulé Accurate characterization of Monte Carlo calculated electron beams for radiotherapy, Medical Physics, 24, pages 401-416, 1997, décrit un accélérateur linéaire médical délivrant un faisceau de photons de haute énergie (quelques Méga-Volts, MV). Trois groupes particuliers de photons se 30 distinguent : les photons provenant de la cible de conversion sur laquelle frappent les électrons incidents, les photons ayant diffusé dans le collimateur 305 1 9 3 7 3 primaire et les photons ayant diffusé dans le cône égalisateur. Plusieurs stratégies ont été proposées pour construire les modèles de sources virtuelles : (a) soit à partir des données contenues dans les espaces des phases pré- 5 calculées par modélisation MC du faisceau d'irradiation, (b) soit à partir de données expérimentales. Le document de A.E. Schach von Wittenau, et al, intitulé Correlated Histogram Representation of Monte Carlo Derived Medical Accelerator Photon-Output Phase Space, Medical Physics, 26 (7), pages 10 1196-1211, 1999, décrit l'utilisation d'histogrammes multidimensionnels corrélés qui décrivent les corrélations entre les variables physiques décrivant l'état du faisceau, à savoir l'énergie, la position et la direction de chaque particule. Le nombre de ces variables est différent suivant le repère choisi et le modèle de source à construire. Il y a un histogramme corrélé par source 15 élémentaire. L'avantage de cette approche est de conserver au sein des histogrammes les corrélations existantes entre toutes les variables stockées dans le fichier d'espace des phases, ce qui permet une représentation plus compacte. La publication de M. Fix, et al, intitulé Monte Carlo Source Model 20 for Megavoltage Photon-Beam Radiotherapy: Photon Source Characteristics, Medical Physics, 31 (11), pages 3106-3121, 2004, divulgue un système dans lequel sont représentés les trois composants physiques essentiels dans la tête de l'accélérateur linéaire médical : la cible, le collimateur primaire, et le cône égalisateur. La population de photons issus de ces trois composants 25 possède des caractéristiques en position, en direction et en énergie bien particulières. Le procédé décrit va chercher à représenter avec le plus de précision possible le comportement physique de ces différentes populations de photons, à différentes énergies (6 et 18 MV). Les auteurs réalisent donc un calcul de l'espace des phases par méthode MC à l'aide d'un code MC 30 appelé BEAMnrc et connu de l'homme du métier dans le domaine de la radiothérapie, pour un accélérateur linéaire médical donné. Le fichier 3051937 4 d'espace des phases est ensuite stocké au-dessus des systèmes de collimation secondaires (mâchoires et collimateur multi-lames). Un premier tri des particules est réalisé pour séparer les photons des électrons, ces derniers étant stockés dans un fichier à part. Les photons constituent la 5 population majoritaire de l'espace des phases et sont eux-mêmes triés en trois groupes, selon que leur dernière interaction a eu lieu dans la cible, le collimateur primaire ou le cône égalisateur. Les auteurs définissent ensuite un repère (cf Figure 1), dans lequel les positions et directions des photons vont être représentées. Dans ce repère, trois plans représentés sur la Figure 10 1 servent à mettre en place le modèle pour chaque source élémentaire. Le premier plan 1, noté « plan du PSF », est celui dans lequel est enregistré l'espace des phases lors du calcul MC préliminaire. Le second plan 2 est celui du plan d'échantillonnage : il est positionné de manière à être confondu avec le plan de l'isocentre, mais peut être a priori localisé ailleurs. 15 L'altitude, par rapport au plan du PSF, du troisième plan 3 dépend de la source élémentaire modélisée, car il correspond au plan de sortie du composant qu'elle représente. Il est appelé « plan origine ». Les coordonnées de chaque photon stockées dans le fichier d'espace des phases, au niveau du plan du PSF, sont notées (x, y). Les coordonnées sont 20 projetées dans le plan origine et le plan d'échantillonnage par « lancer de lumière » ou en anglo-saxon « ray-tracing », en utilisant la direction des photons (u,v), elles sont également stockées dans le fichier d'espace des phases. Les coordonnées polaires de la particule dans le plan origine et dans le plan d'échantillonnage sont notées (ro, cpo) et (rs, cos), respectivement. La 25 différence des angles cps et (po est notée (Prel. L'espace des phases de chacune des trois sources élémentaires est alors décrit par une fonction à plusieurs dimensions qui peut être modélisée par les équations 1 ou 2 (ci-dessous), où fi représente la position et la direction des particules appartenant à la source élémentaire i et gi son 30 énergie. Ces fonctions sont décrites par des histogrammes qui dépendent du modèle choisi. 3051937 5 Fi(x, y , v , E) = fi(rs, (ps,ro(rs), (Po(rs)). (rs)) (équation 1) Fi(x, y , , E) = fi(rs,cps,ro(rs), cf) 0 (rs, ro)). g i(E (rs, ro)) (équation 2) E : énergie 5 A chacune de ces équations correspond en fait un modèle de source. Le premier modèle, correspondant à l'équation 1, est inspiré de celui proposé par Deng et al, Photon Beam Characterization and Modelling for Monte Carlo Treatment Planning, Physics in Medicine and Biology, 45 (2): pages 411427, 2000. Ce modèle de source compact, abrégé MSC par la suite, ne 10 prend en compte que la dépendance de ro, E, et (po avec rs. L'angle (Prel traduit la direction de la particule, c'est donc lui qui est utilisé dans le modèle. On construit tout d'abord un histogramme principal qui correspond à la distribution de la variable rs. Pour chacun des ntot canaux de l'histogramme, on détermine les histogrammes secondaires associés en ro, E, et (Orel. Ces 15 histogrammes représentent les distributions en ro, E, et (Orel des particules p qui ont une position radiale rs, comprise dans le nième canal de l'histogramme principal. On a donc autant d'histogrammes secondaires en ro, E, et (Orel que de canaux dans l'histogramme principal en rs, soit ntot. Le second modèle (équation 2) développé par Fix, propose une amélioration par rapport au 20 premier modèle, qui consiste à prendre en compte une dépendance supplémentaire, celle de l'énergie E et de l'angle (po avec ro. Comme pour le premier modèle, on construit un histogramme principal qui correspond à la distribution de la variable rs. On construit également un histogramme secondaire en ro pour chaque canal de l'histogramme principal. Pour chacun 25 des mtot canaux des ntot histogrammes secondaires en ro, on construit les histogrammes des variables E, et (Orel. Ces histogrammes décrivent respectivement la distribution en énergie et en angle (Orel des particules p dont les positions radiales rs, et ro, sont comprises dans le mième canal de l'histogramme secondaire en ro correspondant au nième canal de 30 l'histogramme principal en rs. On a donc un histogramme en rs de ntot 305 1 9 3 7 6 canaux, ntot histogrammes en ro composés de mtot canaux, et enfin ntot x mtot histogrammes en E, et (prel. L'utilisation des modèles se fait par des tirages aléatoires successifs dans les différents histogrammes. La précision de chacun des deux modèles est évaluée source par 5 source, grâce à la comparaison de distributions de dose dans l'eau obtenues avec les modèles de source et par simulations MC conventionnelles. Plusieurs tailles de champ (1 x 1 cm2 à 30 x 30 cm2) et plusieurs distances source-profondeur dans la cuve à eau (50, 75, 100 et 200 cm) sont testées. La taille des voxels est de 0,4 x 0,4 x 0,4 cm3. Le code utilisé est alors 1 0 DOSXYZnrc connu de l'homme du métier. La simulation MC est ici considérée comme étant la référence. Cette étude est l'une des seules à proposer une validation de MSC pour des distances source-profondeur supérieures à 200 cm, ce qui est intéressant puisqu'il s'agit de la distance possible des imageurs embarqués. Le modèle 1 produit de bons résultats 15 pour une distance de 50 cm, la supériorité du modèle 2 permet de mieux reproduire le faisceau réel que le modèle 1. L'invention concerne un procédé pour estimer les paramètres caractéristiques physiques ou variables d'un faisceau de corpuscules destiné à être émis vers un corps ou un objet et pour déterminer la dose déposée par 20 le faisceau émis, les corpuscules étant rangés dans un fichier d'espace des phases, caractérisé en ce qu'il comporte au moins les étapes suivantes : - Collecter les corpuscules présents dans un faisceau ayant traversé un collecteur et un filtre, ayant interagi ou non avec eux, - Déterminer un nombre d'intervalles I consécutifs et contigus pour le 25 classement des corpuscules du fichier PSF en énergie, position et direction, la réunion de ces intervalles pour chaque variable du faisceau constituant un canal multi-varié, - Découper un intervalle I en un nombre Dmax de micro-canaux réguliers de taille donnée Ô et calculer les vraisemblances correspondant aux 30 intervalles possibles par concaténation contigüe des micro-canaux, retenir le meilleur partitionnement de I, 3051937 7 - Définir un maillage dans chaque direction de l'espace des phases en maximisant la vraisemblance de la réalisation d'une variable dans un intervalle en introduisant un « a priori », afin d'obtenir une partition régularisée de l'espace des phases, 5 - Estimer la distribution de probabilité de l'espace des phases pour la partition régularisée, - Calculer les caractéristiques physiques du faisceau de corpuscules, les comparer à des données prévues et déterminer les doses avant d'émettre le faisceau. 10 L'invention concerne aussi un dispositif pour estimer les paramètres caractéristiques physiques d'un faisceau de corpuscules destiné à être émis vers un corps ou un objet et pour déterminer la dose à émettre, les corpuscules étant rangés dans un fichier d'espace des phases, caractérisé en ce qu'il comporte au moins les éléments suivants : 15 - Un collecteur des corpuscules émis dans un faisceau, - Un processeur adapté à exécuter les étapes du procédé selon l'invention. Le dispositif est par exemple, un dispositif médical et les corpuscules les rayons X pour lesquels on souhaite déterminer la dose à émettre. 20 D'autres caractéristiques et avantages de la présente invention apparaîtront mieux à la lecture de la description d'exemples donnés à titre illustratif et nullement limitatifs annexée des figures qui représentent : - Figure 1, une représentation des plans choisis par l'art antérieur pour la représentation d'histogrammes de distribution de particules, 25 - Figure 2, un exemple d'architecture d'un système accélérateur de particules (dit aussi accélérateur linéaire médical), - Figure 3, une représentation des densités de probabilités de variables et leur densité, - Figure 4, une représentation de l'application du mécanisme d'arbre de 30 Polya. 3051937 8 Le procédé et le système selon l'invention, pour caractériser et déterminer une dose déposée par un faisceau de corpuscules, reposent notamment sur la mise en oeuvre d'un modèle de source compact (MSC) basé sur des histogrammes multi-variés de la distribution des particules dans 5 un faisceau, de dimension supérieure à trois prenant en compte toutes les corrélations possibles entre les variables de l'espace de phases. La taille mémoire à utiliser pour stocker par exemple les caractéristiques du faisceau est réduite selon le procédé de MSC. Tout en réduisant la taille mémoire nécessaire, le modèle de source contient toute l'information présente dans le 1 0 PSF, en particulier les corrélations entre toutes les variables nécessaires à la caractérisation du faisceau (énergie, position, direction des particules). La réduction de la taille mémoire facilite l'intégration logicielle du modèle et accélère sa manipulation par rapport au PSF car le MSC peut être entièrement stocké en une fois en mémoire vive de l'ordinateur. Le procédé 15 selon l'invention consiste notamment à utiliser un maillage irrégulier des canaux, optimal au sens d'un critère statistique de type maximum a posteriori et une technique de lissage des histogrammes produits agissant comme un débruitage de ces derniers (histogrammes qui comportent les paramètres physiques caractéristiques d'un faisceau, à savoir l'énergie, la position et la 20 direction des particules contenues dans un faisceau). D'autre part, la mise en oeuvre du procédé permet d'utiliser un PSF pour la création du MSC moins volumineux que si le calcul de dose était effectué en l'absence de MSC, soit directement à partir de ce dit PSF. Le procédé permet ainsi de réduire le temps de simulation requis pour la constitution du PSF et augmente la 25 précision des calculs de dose. La figure 2 donne un exemple d'architecture d'un LINAC comprenant un module de mise en oeuvre des étapes du procédé selon l'invention. Le système accélérateur de particules (dit aussi accélérateur 30 linéaire médical) 20 comporte une source primaire de particules 21, typiquement des électrons de haute énergie (MV), issus d'un canon à 3051937 9 électrons et accélérés dans un tube où règne le vide (appelé section accélératrice) par l'onde électromagnétique produite par un électro-aimant. Le faisceau d'électrons primaires rencontre ensuite une cible de conversion en matériau lourd 22, produisant ainsi un faisceau divergent et multi-énergies 5 de rayons X à haute énergie par interaction (Bremsstrahlung) des électrons avec la cible. Le faisceau de rayons X est ensuite mis en forme, en termes de géométrie et d'intensité, à l'aide de différents systèmes de collimation, collimateur primaire 23, et d'égalisation, filtre égalisateur 24. On retrouve ensuite d'autres éléments de collimation, par exemple des collimateurs multi- 10 lames 25 et des mâchoires 26, qui permettent de conformer le faisceau d'irradiation pour chaque traitement, individuellement. Les composants 21, 22, 23 et 24 constituent la partie haute de la tête de l'accélérateur, qui reste fixe pour toutes les irradiations. Sur les nouvelles générations d'accélérateurs de particules, le composant 24 tend à disparaître. Les 15 composants 25 et 26 sont mobiles et leur placement, complexe, dépend du traitement à délivrer et est calculé au cas par cas. Le PSF peut alors être calculé en sortie de la partie haute du système accélérateur de particules (sortie du composant 24) en utilisant un processeur 27 comportant un algorithme de traitement adapté à exécuter les étapes du procédé selon 20 l'invention et fournir le MSC. Le processeur 27 comporte également un module de calcul de dose utilisant le MSC. La dose calculée peut alors, par exemple, être utilisée pour optimiser le faisceau de corpuscules émis par le dispositif médical pour le traitement des tumeurs. Les faisceaux de rayons X utilisés en imagerie médicale sont de 25 plus basse énergie (kV) et issus d'un tube à rayons X. Dans un tube à rayons X, une cible en matériau lourd (appelée anode) est bombardée par un faisceau d'électrons accélérés dans le vide, ces électrons étant obtenus en chauffant un filament appelé cathode (effet thermo-ionique) et en étant accélérés par une forte différence de potentiel. Des filtres additionnels 30 (cuivre, aluminium) sont utilisés pour mettre en forme le spectre en énergie des photons issus de la cible. 3051937 10 Le fichier PSF représente un échantillon statistique des particules ou corpuscules représentant la distribution des particules dans le faisceau. Le procédé consistant à déterminer les caractéristiques du faisceau selon un histogramme multidimensionnel (ou multi-variés) consiste à ranger chaque 5 particule présente dans le PSF dans un canal multi-varié. Ainsi, pour chaque variable intervenant dans la caractérisation du faisceau (énergie, position, direction des particules), des intervalles consécutifs et contigus doivent être déterminés pour permettre le classement des particules ou corpuscules. La réunion de ces intervalles pour chaque variable du faisceau constitue un 10 canal multi-varié. Le modèle de source consistant à représenter la distribution des particules du faisceau selon un histogramme multi-varié revient donc à considérer que la densité de probabilité de cette distribution est uniforme sur chaque canal considéré. Le procédé selon l'invention va utiliser deux étapes principales 15 consistant pour la première à réaliser un maillage irrégulier, adaptatif et régularisé de l'espace multidimensionnel dans lequel le faisceau de particules est observé, pour la deuxième à utiliser une technique d'estimation de densité de probabilité qui peut être interprétée comme un lissage des histogrammes corrélés obtenus après accumulation des particules dans les 20 canaux issus du maillage. L'issue de ces deux étapes permettra de déterminer la dose déposée par les particules émises. Le processeur 27 va donc être adapté à exécuter les étapes du procédé détaillées ci-après. L'optimisation de la largeur des intervalles est traitée 25 indépendamment pour chaque variable intervenant dans la source compacte : énergie, rayon ou positions en l'absence de symétrie radiale, angle polaire, angle azimutal. Le maillage multidimensionnel est ainsi séparable en chaque variable intervenant dans le même modèle. Dans ce qui suit, la technique est présentée pour une des variables considérées 30 (énergie, position, directions). 3051937 11 La première étape est de créer des canaux irréguliers pour le rangement des particules du PSF dont les largeurs sont optimisées par un critère du maximum de vraisemblance, par exemple. Ceci est réalisé en appliquant une modélisation probabiliste des données stockées sous forme 5 d'histogramme pour chaque variable précitée intervenant dans la caractérisation du faisceau. La figure 3 schématise la densité de probabilité théorique d'une variable arbitraire de l'espace des phases, 30. L'estimation de cette distribution par le MSC pris dans la dimension de la variable est présentée 10 sous la forme d'histogramme normalisé, 31. La densité du MSC est constante sur un canal. Le maillage représenté dans cet exemple est irrégulier et représentatif du résultat de l'algorithme adaptatif utilisé. La log-vraisemblance, notée L1, est dans le cas d'histogrammes irréguliers contenant les échantillons d'une des variables caractérisant les 15 particules du PSF définie par l'équation (1) : L1 A log(P(XII)) (1) où X est une variable aléatoire décrivant la répartition des particules du PSF au sein de l'ensemble des intervalles I. Cet ensemble est partitionné en une succession de D canaux dont le lième a pour largeur II. Dans le cas d'une source à 4 variables, D peut être typiquement fixé à 64. Si n et N1 désignent 20 respectivement le nombre total de particules totales du faisceau et le nombre de particules dans l'intervalle II, alors X = {A/1, N2, ..., ND}, I = 12, , ID} et P(XII), la probabilité de la réalisation de la variable X sachant la collection d'intervalles I s'écrit : D N (2) 13 (XII) = (N j=1 nIIII Cette expression est caractéristique du modèle de source proposé pour 25 lequel la densité de probabilité des variables caractérisant le faisceau est considérée constante par canal et revient à la probabilité d'avoir N1 particules du PSF dans l'intervalle 305 1 9 3 7 12 L'histogramme empirique associé à la partition en intervalles I présente une Ni , densité de probabilité estimée pour le canal j (de largeur Ii) égale à 7n . Des lors, la probabilité d'observer N1 particules indépendantes venant de ce canal ( N- Ni )1 est = . Les D canaux étant indépendants, on obtient l'expression de nlIil 5 P(Xli), la probabilité d'observer N particules indépendantes venant de l'ensemble des canaux. On en déduit l'expression de la log-vraisemblance des intervalles en fonction des particules du PSF : D Nj ) L1= Nj. log nifil j=1 L'expression de (-L1) peut être vue comme une fonction de coût à 10 minimiser. On utilise par exemple un critère quantitatif qui permet de contrôler l'approximation faite en représentant l'espace des phases par un histogramme. Le maillage retenu est celui qui est optimal au sens de ce critère. L'algorithme présenté ci-dessous s'appuie sur la programmation 15 linéaire pour maximiser le critère L1 donné par exemple dans la publication de Rozenholc et aL, Combining regular and irregular histograms by penalized likelihood, Computational Statistics and Data Analysis, (2010), 54, 33133323. Le processus de maximisation suit les trois étapes : (a) L'intervalle I est découpé en un grand nombre Dm' (typiquement 20 Dm' = 2048) de micro-canaux réguliers de tailles égales à ô. Les vraisemblances correspondant à tous les intervalles possibles constitués par la concaténation contigüe des micro-canaux sont calculées ; (b) Ces micro-canaux sont associés selon tous les regroupements contigus possibles, à chaque fois de manière à produire D canaux ; (3) 3051937 13 (c) La meilleure association, et donc le meilleur partitionnement de I, est celui qui maximise L1. C'est ce maillage qui sera retenu dans le procédé selon l'invention. Une implémentation possible de ces trois étapes est présentée ci- 5 dessous Étape (a) : Pour i = 0 ... Dmax Dmax : nombre de micro-canaux donnés initialement Initialiser M = 0 Pour j = i + 1 ... Dm' boucle sur les micro-canaux j supérieurs à i. Calculer M = M + rnj où m1 représente le nombre de particules dans le micro-canal 1. Calculer pi(i, j) =1'Lin. log num_06 probabilité de regrouper les micro-canaux de i à j pour former un canal. Étape (b) : Initialiser pour j = 1 ... Dmax boucle sur les j micro-canaux qU) = Pour i = 2 ... D (nombre de canaux du PSF) Pour j = i ...Dm' boucle sur les j micro-canaux Initialiser S = Pour c/ = - 1...j - 1 s = q(d) + pi(d,j) Si s > S S = s H(i,j) = p(i, j) = q(H(i,j)) + p1(H(i,]),]) Pour j = i ...Dm' boucle sur les j micro-canaux q(j) = p(i, j) Étape (c) : Initialiser L = Dmax avec Dmax : nombre de micro-canaux donnés initialement Pour j = D ...1 10 15 20 25 30 3051937 14 L = c(H(j , L[1]), L) où L[1] représente la première coordonnée du vecteur L et c(u, L) la concaténation des vecteurs u et L, on obtient le partitionnement des intervalles. 5 L est ainsi le vecteur contenant les indices de micro-canaux correspondant aux frontières de la partition optimale au sens ML constituée de D canaux. L'étape suivante consiste à générer un maillage régularisé et adaptatif afin de gérer plus facilement les histogrammes corrélés de grandes dimensions et donc de gagner en temps de calcul. Ce maillage va 10 notamment permettre de décrire aussi précisément que possible les détails des distributions en énergie, en position et en direction, déterminant ainsi la précision de la source compacte. Maillage adaptatif régularisé. Cette méthode de maillage est inspirée de celle proposée par 15 Rozenholc et al. Elle sera désignée par l'acronyme « MAP » pour Maximum a posteriori. L'intérêt est de limiter les grandes disparités entre les largeurs des canaux qui peuvent apparaître quand le critère ML seul est utilisé. Dans le cas des fichiers de l'espace des phases PSF faiblement à modérément peuplés, le critère ML est alors sensible à la statistique limitée et peut 20 présenter une tendance à produire un maillage bruité. Cette tendance tend à disparaître avec des PSF fortement peuplés. Afin de proposer un maillage tolérant à une statistique de PSF limitée, le procédé comporte une étape qui régularise les solutions. Le principe de cette méthode consiste à maximiser à nouveau la vraisemblance, mais cette fois-ci, après introduction d'un a priori 25 qui favorise des aires de canaux égales, i.e., une répartition en quantiles. L'a priori choisi prend ici la forme d'une loi de Dirichlet symétrique par exemple. D'après le théorème de Bayes, en conservant les notations de l'approche ML, il est possible d'écrire la probabilité de réalisation de la variable X dans l'intervalle I de la façon suivante : P (X , I) = P(IIX)P (X) = P (XII)P (I) (4) 30 De cela découle l'expression de la vraisemblance pénalisée L2 : 3051937 15 L2 A log(P(XII).P(I)) (5) où P(I) désigne l'a priori modélisé par une loi de Dirichlet symétrique. Sa densité de probabilité s'exprime de la façon suivante : P(I) r T(D a) rif-1 N2 a-1 ND a-1 (6) (a)D n n n où a est un paramètre positif ajustable par l'utilisateur et r(x) représente la fonction Gamma. 5 L'expression de la vraisemblance pénalisée est finalement donnée par l'équation (7) : D D L2 = Ni. log (-n.(1.1) + (a -1).log (A:7j j=1 i=i Le choix de valeur de a permet de pondérer l'importance que l'on donne à l'a priori. Il est ainsi possible de favoriser un maillage proche des quantiles ou la méthode ML. Si a = 1, le maillage choisi est celui de la méthode ML. Si a est 10 très grand, on se rapproche de celui des quantiles. Typiquement a est choisi égal à 100 par des considérations heuristiques. En effet, l'expression de L2 donnée par l'équation (7) suggère que le paramètre a intervient en substance comme un nombre de particules vis-à-vis de la première somme. Ainsi, il est possible de considérer en première approximation que pour des 15 canaux peuplés avec moins de a (100) particules, l'influence de l'a priori est significative. Le maximum de la vraisemblance pénalisée est obtenu en suivant les mêmes étapes que celles décrites pour l'approche ML. L'unique différence réside dans le calcul de pl(i,j) à l'étape (a). 20 Pour le maillage régularisé on calcule pour tout i,j : = log nu_06+ (a - 1) log 11n. Une fois le maillage dans chaque dimension réalisé, une partition régularisée de l'espace des phases est obtenue. L'étape suivante consiste à estimer la distribution de probabilité de l'espace des phases pour cette 25 partition fixée. (7) 3051937 16 L'usage est, par exemple, de comptabiliser dans un canal particulier de l'histogramme multidimensionnel les particules dont les valeurs de variables associées sont comprises dans les intervalles correspondants au canal. Cette façon de comptabiliser des particules est connue en 5 statistiques mathématiques comme l'obtention d'un histogramme empirique. L'histogramme empirique après normalisation, correspond à l'estimateur par maximum de vraisemblance de la distribution de probabilité de l'espace des phases pour la partition de l'espace donnée (cf. Rozenholc et al. précité). La taille de l'échantillon de particules intervient également ici comme élément 10 limitant de cet estimateur empirique. En effet, chaque canal est alors peuplé suivant la réalisation d'une loi de Poisson dont l'intensité est croissante avec le nombre de particules dans le PSF. Or, plus cette intensité est faible plus le rapport signal/bruit (Moyenne Poisson m = À ; Variance Poisson a2 = 1; d'où B -S = -m = .) du canal en question l'est également. Si les histogrammes 1D a 15 correspondant aux lois marginales de chaque grandeur sont largement peuplés, il n'en est pas du tout de même pour une loi jointe en quatre ou cinq dimensions. En effet, dans ce cas, même un espace de phases présentant plus de 109 particules demeure significativement creux, ce qui se traduit par des canaux d'histogrammes multidimensionnels très bruités. Cet effet est 20 une conséquence directe de la grande dimensionnalité associée à la modélisation physique d'un système complexe, tel un accélérateur ou un tube à rayons X. Il devient donc légitime de considérer, à l'image de la démarche visant à régulariser la taille des canaux rappelée plus haut, la possibilité d'injecter une 25 connaissance a priori sur la forme des distributions d'espaces de phases afin de modérer le bruit induit par la taille limitée du PSF. Le procédé comporte une étape complémentaire consistant à lisser les histogrammes obtenus afin de limiter les artefacts liées à une statistique limitée se traduisant par des fluctuations aléatoires induites dans 30 les canaux (bruit Poisson). L'étape de lissage du procédé consiste en la réduction de ces fluctuations statistiques dans les canaux du MSC. 3051937 17 La figure 4 illustre le fonctionnement d'un estimateur par lissage par arbre de Pôlya. L'arbre de Pôlya est une distribution de probabilité. Pour illustrer de manière simple l'étape de lissage mise en oeuvre par le procédé selon l'invention, l'exemple qui suit est donné dans un premier temps, pour le 5 cas d'un arbre de Pôlya monodimensionneL Il est possible d'imaginer un arbre binaire partitionnant l'espace et une collection de nombres positifs a associés à chaque branche de l'arbre. La distribution de l'arbre de Pôlya est obtenue en générant récursivement des variables aléatoires de loi Beta paramétrée par les nombres a affectés à 10 chaque sous-branche correspondant à la ramification de la branche de niveau supérieur. La masse de probabilité associée à une branche est répartie aléatoirement sur deux sous-branches. Les notations doivent être choisies afin de pouvoir identifier chaque niveau de ramification de l'arbre. Pour ce faire, l'indice du canal de la partition correspondant à une branche 15 est représenté par un codage binaire de celui-ci, ce qui permet d'identifier clairement chaque niveau. Ainsi, pour un niveau m particulier, on note a''...,', avec ci = 0 ou ci = 1 le paramètre associé au canal (branche) k noté E1E2 ...Em en codage binaire (à base 2). La ramification de cette branche (découpe aléatoire du canal B£1£2 ...£m en deux sous-parties /3,1,2...Emo et 20 BE1E2...Emi), est alors obtenue par la génération d'une variable de loi Beta(aE1E2...Emo, aE1E2...,m1). La procédure de ramification de l'arbre peut être poursuivie indéfiniment. En pratique, on l'arrête bien entendu à un niveau prédéfini M. Il est fréquent de choisir des coefficients a£1£2 ...£m égaux pour un niveau m donné. En particulier, un choix usuel revient à prendre aEiE2 ...Em = 25 am = c m2 ou am = c 2m, où c est un nombre positif (typiquement c = 1). Cette construction récursive est répétée jusqu'au niveau de ramification maximal souhaité M. Cet a priori présente un certain nombre d'avantages. Notamment, l'expression de la loi a posteriori de l'arbre de Pôlya, prenant en compte les 30 données de l'espace des phases, peut être obtenue de façon analytique. Cet 3051937 18 élément est déterminant dans la mesure où une source compacte (en 4D ou plus) fait émerger un nombre de canaux très important, ce qui empêche en pratique de concevoir un algorithme de maximisation en vue d'optimiser un critère statistique de type maximum a posteriori (a minima de l'ordre de 107 5 paramètres). De ce fait, l'expression analytique de l'espérance a posteriori de la distribution à estimer qui, dans le cas monodimensionnel, peut s'écrire pour le canal /3,1,2...Ein : + nEl a + n FT g)(x E = ao +a1+ n + + est le nombre de particules de l'espace des phases présentes dans le canal /3,1,2...,1 et n le nombre total de particules du PSF. On peut 10 remarquer que si tous les coefficients a sont nuls, on trouve 3)(x E = nE En2 Em , c'est-à-dire l'estimateur ML (empirique) de la source, l'a priori n'intervient alors pas. Dans le cas de sources multidimensionnelles, le passage à un arbre multidimensionnel requiert un certain nombre d'extensions. La 15 modification majeure réside dans le remplacement de la loi Beta par une loi de Dirichlet (voir Annexe) permettant une découpe en plusieurs dimensions [Jara, A.; Hanson, T. E. & Lesaffre, E. Robustifying Generalized Linear Mixed Models Using a New Class of Mixtures of Multivariate Polya Trees, Journal of Computational and Graphical Statistics, 2009, 18, 838-860.]. Comme il a été 20 indiqué, de par sa forte dimensionnalité un modèle histogrammé de source compacte retenant toutes les dimensions de l'espace des phases (4D dans le cas d'une symétrie axiale, 5D sinon) présente un nombre de canaux très important. Par exemple, si on limite en pratique le nombre de canaux à 64 pour chaque variable (M = 6 niveaux de ramification par dimension), le 25 nombre de canaux total du modèle de source est alors égal à 644 = 16777216. Dans la suite, l'estimateur de la distribution de l'espace de phases est exprimé pour une dimension arbitraire K et une profondeur d'arbre M. m J 3051937 19 Afin de maintenir des expressions relativement compactes, les notations (k) (k suivantes sont posées : Ek) tm = (k)E2( Em) pour tout 1 REVENDICATIONS1 - Procédé pour estimer les paramètres caractéristiques physiques ou variables d'un faisceau de corpuscules destiné à être émis vers un corps ou un objet et pour déterminer la dose déposée par le faisceau émis, les corpuscules étant rangés dans un fichier d'espace des phases, caractérisé en ce qu'il comporte au moins les étapes suivantes : - Collecter les corpuscules présents dans un faisceau ayant traversé un collecteur (23) et un filtre (24), ayant interagi ou non avec eux, - Déterminer un nombre d'intervalles I consécutifs et contigus pour le classement des corpuscules du fichier PSF en énergie, position et direction, la réunion de ces intervalles pour chaque variable du faisceau constituant un canal multi-varié, - Découper un intervalle I en un nombre Dmax de micro-canaux réguliers de taille donnée Ô et calculer les vraisemblances correspondant aux intervalles possibles par concaténation contigüe des micro-canaux, retenir le meilleur partitionnement de I, - Définir un maillage dans chaque direction de l'espace des phases en maximisant la vraisemblance de la réalisation d'une variable dans un intervalle en introduisant un « a priori », afin d'obtenir une partition régularisée de l'espace des phases, - Estimer la distribution de probabilité de l'espace des phases pour la partition régularisée, - Calculer les caractéristiques physiques du faisceau de corpuscules, les comparer à des données prévues et déterminer les doses avant d'émettre le faisceau. 2 - Procédé selon la revendication 1 caractérisé en ce que pour définir le maillage on exécute les étapes suivantes : (a) découper l'intervalle I en un grand nombre Dmax de micro-canaux réguliers de tailles égales à 8, 3051937 24 (b) associer ces micro-canaux selon tous les regroupements contigus possibles à chaque fois de manière à produire D canaux, (c) déterminer le partitionnement de I, qui maximise la valeur L1, log vraisemblance des intervalles en fonction des corpuscules du fichier PSF 5 L1 = E13_,Ni.log(), où D est le nombre de canaux du fichier PSF, Nj le nli j I nombre de corpuscules observables du canal j, lj l'intervalle pour le canal j, n le nombre total de corpuscules dans le faisceau. 3 - Procédé selon la revendication 2 caractérisé en ce qu'il comporte pour 10 optimiser la valeur L1 au moins les étapes suivantes : Étape (a) : Pour i = 0 ... Dmax Dmax : nombre de micro-canaux donnés initialement Initialiser M = 0 Pour j = i + 1 ...Dm' boucle sur les micro-canaux j, j>i, Calculer M = M + rnj où m1 représente le nombre de particules dans le micro-canal 1. Calculer pi(i, j) =11n. log num_06 probabilité de regrouper les micro-canaux de i à j pour former un canal, Étape (b) : Initialiser pour j = 1 ...Dm' boucle sur les j micro-canaux qU) = Pour i = 2...D (nombre de canaux du PSF) Pour j = i ...Dm' boucle sur les j micro-canaux Initialiser S = Pour c/ = - 1...j - 1 s = q(d) + pi(d,j) Si s > S S = s H(i,j) = 15 20 25 3051937 25 p(i,j) = q(H(i, j)) + pi(H(i, j), j) Pour j = i ... Dm' boucle sur les j micro-canaux a(j) = p(i,j) Étape (c) : Initialiser L = Dmax Dmax : nombre de micro-canaux donnés initialement Pour j = D ...1 L = c(H(j , LW), L) où L[1] représente la première coordonnée du vecteur L et c(u, L) la concaténation des vecteurs u et L, on obtient le partitionnement des intervalles. 4 - Procédé selon l'une des revendications 1 à 3 caractérisé en ce que l'on utilise une loi de Dirichlet symétrique comme « a priori ». 15 5 - Procédé selon l'une des revendications 1 à 4 caractérisé en ce que pour estimer la distribution de probabilité des corpuscules de l'espace des phases PSF pour une partition régularisée, on comptabilise dans un canal particulier de l'histogramme multidimensionnel formé par intervalles pour les différentes variables, énergie, position, direction des corpuscules, les corpuscules dont 20 les valeurs de variables associées sont comprises dans les intervalles correspondants. 6 - Procédé selon la revendication 5 caractérisé en ce qu'il comporte une étape de lissage des histogrammes obtenus en utilisant un arbre de Pôlya 25 multidimensionnel et en exprimant la loi a posteriori de l'arbre de Pôlya selon une expression analytique. 7 - Procédé selon l'une des revendications précédentes caractérisé en ce que l'on détermine une dose déposée par des rayons X d'un dispositif 30 médical. 5 10 3051937 26 8 - Dispositif pour estimer les paramètres caractéristiques physiques d'un faisceau de corpuscules destiné à être émis vers un corps ou un objet et pour déterminer la dose à émettre, les corpuscules étant rangés dans un fichier d'espace des phases, caractérisé en ce qu'il comporte au moins les 5 éléments suivants : - Un collecteur (23) des corpuscules émis dans un faisceau, - Un processeur (27) adapté à exécuter les étapes du procédé selon l'une des revendications 1 à 7. 10 9 - Dispositif selon la revendication 8 caractérisé en ce que le dispositif est un dispositif médical et les corpuscules des rayons X.