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 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 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,

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

L'article de CM. Ma, B.A. Faddegon, D.W.O. Rogers et T.R.

Mackie intitulé Accurate characterization of Monte Carlo calculated électron beams for radiotherapy, Médical Physics, 24, pages 401 -41 6, 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 distinguent : les photons provenant de la cible de conversion sur laquelle frappent les électrons incidents, les photons ayant diffusé dans le collimateur 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é- 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 Représentation of Monte Carlo Derived Médical Accelerator Photon-Output Phase Space, Médical Physics, 26 (7), pages 1 196-121 1 , 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 é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 for Megavoltage Photon-Beam Radiotherapy: Photon Source Characteristics, Médical Physics, 31 (1 1 ), 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 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 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 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 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 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. 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 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 ( ,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 (r _o , φ _o ) et (r _s , φ _s ), respectivement. La différence des angles φ _s et φ _o est notée φ _rel .

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 énergie. Ces fonctions sont décrites par des histogrammes qui dépendent du modèle choisi.

E : énergie

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 41 1 - 427, 2000. Ce modèle de source compact, abrégé MSC par la suite, ne prend en compte que la dépendance de r ₀ , E, et φ _o avec r _s . L'angle φ _rel 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 r _s . Pour chacun des n _m canaux de l'histogramme, on détermine les histogrammes secondaires associés en r ₀ , E, et φ _rel . Ces histogrammes représentent les distributions en r ₀ , E, et φ _rel des particules p qui ont une position radiale r _s , comprise dans le n' ^eme canal de l'histogramme principal. On a donc autant d'histogrammes secondaires en r ₀ , E, et φ _rel que de canaux dans l'histogramme principal en r _s , soit n _to t- Le second modèle (équation 2) développé par Fix, propose une amélioration par rapport au premier modèle, qui consiste à prendre en compte une dépendance supplémentaire, celle de l'énergie E et de l'angle φ _o avec r ₀ . Comme pour le premier modèle, on construit un histogramme principal qui correspond à la distribution de la variable r _s . On construit également un histogramme secondaire en r _o pour chaque canal de l'histogramme principal. Pour chacun des mtot canaux des n _tot histogrammes secondaires en r ₀ , on construit les histogrammes des variables E, et φ _rel . Ces histogrammes décrivent respectivement la distribution en énergie et en angle φ _rel des particules p dont les positions radiales r _s , et r ₀ , sont comprises dans le m' ^eme canal de l'histogramme secondaire en r _o correspondant au n' ^eme canal de l'histogramme principal en r _s . On a donc un histogramme en r _s de n _tot canaux, n _tot histogrammes en r _o composés de m _tot canaux, et enfin n _tot x m _tot histogrammes en E, et φ _rel . 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 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 cm ² à 30 x 30 cm ² ) 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 cm ³ . Le code utilisé est alors 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 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 .

Le brevet US 6,535,837 décrit une méthode pour représenter de façon compacte l'espace des phases afin de caractériser un faisceau de corpuscules pour le calcul de la dose délivrée au patient dans un système de radiothérapie.

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 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 une partition du support Ω de chaque variable sous forme d'un ensemble / d'intervalles 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 ledit support Ω en un nombre D _max de micro-canaux réguliers de taille donnée δ et calculer les vraisemblances correspondant aux intervalles possibles par concaténation contigue des micro-canaux, retenir la meilleure partition / de Ω,

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

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 :

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

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, • 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 Polya.

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 œuvre d'un modèle de source compact (MSC) basé sur des histogrammes multi-variés de la distribution des particules dans 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 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é 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 direction des particules contenues dans un faisceau). D'autre part, la mise en œuvre 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 précision des calculs de dose. La figure 2 donne un exemple d'architecture d'un LINAC comprenant un module de mise en œuvre des étapes du procédé selon l'invention.

Le système accélérateur de particules (dit aussi accélérateur linéaire médical) 20 comporte une source primaire de particules 21 , typiquement des électrons de haute énergie (MV), issus d'un canon à é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 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- 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 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 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 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 (cuivre, aluminium) sont utilisés pour mettre en forme le spectre en énergie des photons issus de la cible.

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 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 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 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 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 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 (énergie, position, directions).

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 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 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 L _{l t} est dans le cas d'histogrammes irréguliers contenant les échantillons d'une des variables caractérisant les particules du PSF définie par l'équation (1 ) :

où X est une variable aléatoire décrivant la répartition des particules du PSF au sein de l'ensemble des intervalles / définissant une partition du support Ω de la variable. Ce support est partitionné en une succession de D canaux dont le j ^ième a pour largeur I _j . Dans le cas d'une source à 4 variables, D peut être typiquement fixé à 64. Si n et N _j désignent respectivement le nombre total de particules totales du faisceau et le nombre de particules dans l'intervalle alors la

probabilité de la réalisation de la variable X sachant la collection d'intervalles / s'écrit :

Cette expression est caractéristique du modèle de source proposé pour lequel la densité de probabilité des variables caractérisant le faisceau est considérée constante par canal et revient à la probabilité d'avoir Nj particules du PSF dans l'intervalle I _j .

L'histogramme empirique associé à la partition en intervalles / présente une densité de probabilité estimée pour le canal ; (de largeur /,) égale a Des lors, la robabilité d'observer N, particules indépendantes venant de ce canal est Les D canaux étant indépendants, on obtient l'expression de

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 :

L'expression de (-L ₁ ) peut être vue comme une fonction de coût à 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 linéaire pour maximiser le critère L ₁ donné par exemple dans la publication de Rozenholc et al., Combining regular and irregular histograms by penalized likelihood, Computational Statistics and Data Analysis, (201 0), 54, 331 3- 3323. Le processus de maximisation suit les trois étapes :

(a) Le support Ω est découpé en un grand nombre D _max (typiquement Dmax = 2048) de micro-canaux réguliers de tailles égales à δ. Les vraisemblances correspondant à tous les intervalles possibles constitués par la concaténation contigue 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 ;

(c) La meilleure association, et donc la meilleure partition / de Ω, est celle qui maximise L ₁ . 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- dessous

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 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 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 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 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, 1) = P(1\X)P(X) = P(X\1)P(1) (4) De cela découle l'expression de la vraisemblance pénalisée L ₂ :

où P(/) 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 :

où a est un paramètre positif ajustable par l'utilisateur et Γ(χ) représente la fonction Gamma.

L'expression de la vraisemblance pénalisée est finalement donnée par l'équation (7) :

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 α = 1, le maillage choisi est celui de la méthode ML. Si a est 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 L ₂ 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 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 P ₁ (i,j) à l'étape (a). Pour le maillage régularisé on calcule pour tout i,j :

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 partition fixée.

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 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 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 σ ² = λ; d'où du canal en question l'est également. Si les histogrammes 1 D

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 10 ⁹ particules demeure significativement creux, ce qui se traduit par des canaux d'histogrammes multidimensionnels très bruités. Cet effet est 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 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 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.

La figure 4 illustre le fonctionnement d'un estimateur par lissage par arbre de Polya. L'arbre de Polya est une distribution de probabilité. Pour illustrer de manière simple l'étape de lissage mise en œuvre par le procédé selon l'invention, l'exemple qui suit est donné dans un premier temps, pour le cas d'un arbre de Polya 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 Polya est obtenue en générant récursivement des variables aléatoires de loi Beta paramétrée par les nombres a affectés à 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 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 le paramètre associé au canal (branche) k

noté en codage binaire (à base 2). La ramification de cette branche

(découpe aléatoire du canal en deux sous-parties et

est alors obtenue par la génération d'une variable de loi

Βeta 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 égaux pour un niveau m donné. En particulier, un choix usuel revient à prendre 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 données de l'espace des phases, peut être obtenue de façon analytique. Cet é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 10 ⁷ 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

où est le nombre de particules de l'espace des phases présentes

dans le canal et n le nombre total de particules du PSF. On peut

remarquer que si tous les coefficients a sont nuls, on trouve

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 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é 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 nombre de canaux total du modèle de source est alors égal à 64 ⁴ = 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. Afin de maintenir des expressions relativement compactes, les notations suivantes sont posées : pour tout et

représente l'indice d'un canal de la source dans la

dimension k. Par conséquent, un canal d'une source à dimensions sera identifié par les K indices et sera représenté par

Les variables binaires

sont également introduites. Il est ainsi possible de construire un indice . À titre d'exemple, l'expression de l'estimateur d'une

source monodimensionnelle (1 ) devient avec ces notations :

L'espérance a posteriori de l'arbre de Pôlya pour une source multidimensionnelle de dimension K peut alors s'exprimer :

où est le nombre de particules de l'espace des phases

présentes dans le canal Cette expression autorise un calcul rapide de l'estimateur bayésien de source compacte à partir de l'histogramme 4D des particules de l'espace des phases. Ceci peut être interprété comme un post filtrage non linéaire des données histogrammées.

Un choix de paramètres tel que :

permet de centrer l'arbre sur la loi correspondant au produit des lois marginales une dimension, 1 D, des K variables. Ce paramétrage se révèle intéressant notamment quand les différentes variables de la source sont faiblement corrélées. L'arbre de Pôlya n'intervient alors que comme une correction multidimensionnelle conjointe par rapport à une situation où les variables de l'espace de phases seraient considérées indépendantes (à la façon d'une fonction copule).

Une autre solution est dite Arbre de Pôlya « adouci ». L'estimateur de la distribution d'espace de phases présenté ci-dessus est continu en tout point sauf aux frontières de la partition où il peut présenter des discontinuités importantes. Pour pallier ce phénomène, il est possible d'imaginer un modèle de source compacte où la partition n'est plus considérée comme unique mais distribuée de façon à tempérer ces effets de discontinuités. Une manière de procéder consiste à décaler l'origine de la partition multidimensionnelle suivant tous les décalages ou offsets possibles (soit M ^K offsets possibles en dimension K) et de procéder au calcul de l'espérance a posteriori de l'arbre pour chaque décalage. Afin de conserver un a priori identique garantissant la même distribution centrale quel que soit le décalage, les coefficients doivent être modifiés de façon ad hoc. Afin d'expliciter cette

recombinaison nous notons le décalage de la partition pour la dimension

k au niveau M le plus ramifié de l'arbre. Dans la partition décalée, l'indice en codage binaire du canal initial devient

Cette rotation de l'indice suppose implicitement que la distribution marginale de chaque variable présente des valeurs comparables pour les bornes minimum et maximum de son support. Ceci est vérifié intrinsèquement pour les angles et observé pour les variables rayon et énergie.

Pour chaque décalage en K dimensions on reconstruit

un histogramme empirique pour les canaux de la partition

décalée où le nombre de particules de l'espace des phases est alors noté

De la même façon, représente le nombre de particules de

la marginale k dans la partition étendue. En remontant l'arbre jusqu'à sa racine, on peut évaluer pour tout m. Dès lors, les

coefficients de l'arbre permettant de garantir un a priori identique quel que soit Δ s'écrivent :

Il suffit alors de substituer

dans la formule de l'espérance a posteriori de la distribution multidimensionnelle d'espace des phases pour obtenir l'expression de représente le point de l'espace des phases

décalé de A. On note

distribution ainsi obtenue après suppression de l'offset Δ (recalage).

L'estimateur « adouci » est obtenu en prenant la moyenne des espérances a posteriori pour tous les décalages Δ considérés. Avec cette construction, il est possible de prendre en compte M ^K valeurs distinctes de décalages, ce qui n'est pas réaliste d'un point de vue pratique pour une source 4D. On se restreint alors à quelques dizaines de décalages.

Applications étudiées en radiothérapie et imagerie :

- Modèles de sources pour des accélérateurs linéaires médicaux, intégrables aux logiciels utilisant ces modèles (treatment planning

System, TPS),

- Modèles de sources pour des tubes à rayons X équipant des systèmes d'imagerie RX 2D ou 3D,

- Modèles prédictifs de type MC pour les imageurs planaires de type flat-panel (EPID, systèmes d'imagerie cone-beam CT),

- En radiothérapie : dosimétrie hors-champ, aux organes à risque et aux tissus sains recevant une faible dose.

Le procédé et le système selon l'invention permettent notamment de réduire de l'ordre de trois décades la quantité d'information nécessaire à stocker tout en restituant des caractéristiques de faisceau moins sensibles au bruit d'origine statistique.